用于弹性波方程模拟的基于逐元技术的谱元法.pdf

"#+*

!

月

"学术论文"

用于弹性波方程模拟的基于逐元技术的谱元法"

$

林伟军"!王秀明"!!张海澜"

":中国科学院声学研究所!北京"##*#*$:(DSQUO;2AC1;9.!)QQ(!OULC5""!#!

!!B),T;=G3C1C/AEL;321;"#$!)9?2A/10/J4JO

摘要!!由于结合了谱方法和有限元法的优势!谱元法为弹性波方程的数值模拟提供了一种新的有效工具:文中从弹性动力学方程的弱形式出发!详细阐述了使用M;;3<A;正交多项式展开的谱元法4基本理论及相应数学公式:单元内的弹性波场用高阶M/A/3;插值近似!而每个单元上的数值积分44采用P/9??’MCX/22C’M;;3<A;积分方法:如此可得到对角形式的全局质量阵!简化了运算:逐元技术4被用于M;;3<A;谱元程序中!极大地减少了内存和计算需求:最后!四个数值算例被用于验证这种4谱元法的高精度和强适应性:同时证明这种新的谱元法很适合用于那些具有复杂结构!包括起伏自由表面的非均匀介质中的弹性波模拟:

关键词!!谱元!699;<(9多项式!弹性波方程模拟!逐元技术:!!如何利用数值计算的方法更精确地模拟弹性

波的传播一直是国内外研究者的工作重点:随着计算机技术的发展!一些原来影响数值计算方法应用的瓶颈一一被克服!但那些大型的复杂二维或三维问题!对原有的数值模拟方法而言仍然具有挑战性:

当前在弹性波方程模拟中普遍使用的数值方法!主要有有限差分法!有限元法和谱方法:它们有各自的优缺点并在不同的领域得到应用:

有限差分法"可以被看作是一种配置点方N@R$法"它要求近似解在若干个配置点上的余量为零$!但N网格@R没有使用近似解!而直接在配置点"点$上将微分算子用截断的T/1CA级数展开:它的J误差由配置点的数目和T/1CA级数的截断误差决J("

定’:有限差分法由于编程简捷和计算速度较快被广泛应用于计算地球物理和工程地震中:但对于低阶有限差分法!它需要大量的网格点来保证精

##&’#"’$*收稿!$##&’#!’"*收修改稿!$

批准号#"$)中国科学院访问学者基金和联合国教科文组织理论物理基金赞助#"!+#$#!"国家自然科学基金重点项目"

#’./011036/01:0C/:/=:=3!-!.8

($!!

度’:而对于高阶有限差分法!又很难高效地处理(+自由界面和复杂结构问题’:

&(

有限元方法"是一种加权余量法’:在N-R$

加权余量法中!引进了一组可展开的函数"尝试函数或形函数$作为基函数用于对偏微分方程的解作截断展开:为了保证由有限项基函数截断展开定义的近似函数能最佳地逼近偏微分方程的精确解!测试函数被引进用于验证近似解带来的余量是否达到最小:在有限元方法中通常采用P/1;AE03加权余量法!它将测试函数取和尝试函数相同的基!同时将整个求解域分解为若干个子区域"单元$!因此尝试函数被定义在每个单元上!并通过使近

,(

似解在全域上达到最小余量而求得解’:有限元法

适合处理边界和内部复杂结构!在结构分析和瞬

%!*(态模拟方面获得了广泛的应用’!但它在大型的

弹性波方程模拟领域一直没有得到很好应用:其原

>(

!二是因一是有限元法需要耗用大量计算资源’

!

"#

+>

"#(

!而通常的高阶有限低阶有限元存在频散现象’""(

元存在伪波’:

本文采用M;;3<A;正交多项式!由此得到对角4形式的全局质量矩阵:同时为了减少谱元法对内存

"*(,!

和计算量的巨大需求!逐元技术’被首次引入基

谱方法!它将尝试函数取为无限可微的全局函数"如N!将问题从物理空间变换到C9A0;A级数$谱空间进行处理!具有良好的收敛性:理论上谱方法可以在一个最短波长中取最少的网格点而达到

"$(

所需要的精度’:为了处理通常的边界条件!在

于M;;3<A;展开的谱元法中!并分析了由此带来的4内存和计算量的减少:最后!计算了+个实例模型!证明我们这种采用逐元技术的M;;3<A;展开谱元法4很适合用于复杂介质中弹性波的传播模拟:

空间展开时!正交多项式"(G;X=G;ZCAM;;3<A;J4

多项式$经常被用来取代NC9A0;A级数’"

!(:但由于正交多项式的配置点分布很不均匀!给时间步长的选择带来很大的限制

’"+(

:同时!如同有限差分

法一样!谱方法不能自如地处理弯曲的自由界面

和复杂的结构问题:很多方法如弯曲坐标’"+(

和区域分解’

"&(

被用于解决这类问题!但代价是计算量的增加:

谱元法"DO-R$!最早由O/2;A/在">*+年提

出’

",(

!并主要应用于流体动力学:它把有限元法和谱方法相结合!因此兼具了有限元的处理边界和结构的灵活性和谱方法的快速收敛特性:和有限元相比!在达到相同的精度前提下!通过采用较稀疏的单元划分可以减少计算时间和内存需求:D;A0/30等在">>"年首次将谱元法引入到声波方程的数值模拟中!各国研究者在其后做了大

量的相关研究’""!"%+$#(

:但在处理大型计算问题

时!谱元法仍然受到内存和计算量的限制!需要采用特定的技术或在大型计算机上通过并行计算解决

’"*!$"(

:

谱元法的基本思想是在每一个单元上使用谱

方法:同时选取以截断的正交多项式表示的基函数!在各个单元上通过利用配置点插值!以提高级数表示的解的收敛速度:其主要步骤是#""$首先把计算的区域分成许多子域"单元$!每个子域由若干节点"配置点$组成:"$$在每个子域中把近似解表示成截断的正交多项式展开:"!$用P/1;A’03方法求解正交问题的变分格式!得到全局的近

似解:

有两种类型正交多项式被用于谱元法中#分别

为(G;XJ?G;Z正交多项式’"%!"*!$#

(或M;4

;3<A;正交多项式’""!$"

(:它们都是奇异性D29A.’M0C9Z011;问题

的特征函数:但它们的不同特性也导致了两者在谱元法中的不同应用方式:

!!谱元法理论

!%!!弹性波方程及其弱形式

假定"’8%3为我们考察的区域!其中%3为物理空间维数!这样在弹性非均匀介质中!考虑时间(’

#!F(!弹性波方程可以如下表述##

!5!!!L!5M!MN"5!!03"O’#!F(!""

$!5L#5!03$<

5O’#!F(!"$$!5M%ML$

5!03$(5O’#!F(!"!

$其中!5为位移向量!!5M为应力张量!两者关系由广义‘CCE定律定义#!5ML%5MJ"!J!"!这里%5M

J"是弹性张量;$<5和$(5分别为位移和应力边界!#5和$5分别为预知的边界位移和边界应力:加上相应的初始条件!就构成了弹性波方程的强形式表示:

为了用数值的方法解波动方程!运用加权余量原理!得到弹性波方程的弱形式表示#

"&!#

’$N6"&!!$L"&!"$N"&!$$%!"+

$其中!"!$(&!为尝试函数!&!代表尝试函数空间!其成员满足位移边界条件!并在区域内保证一阶导数平方可积:&(’为测试函数!’为测试函数空间!其成员函数在位移边界上值为零!并在区域内保证一阶导数平方可积:"+$式中的内积可如下表示#

"&!#’$L$<!$"&!#!$L$

<!

$#&5!5

$"<"!"&6"&!!$L)

&"5!M$%5J"*"J!"$<"!""

M

,$"&!"$L)&5"5

<

"!

""

%$%+3

"&!$$%

L*)"&$<

%5L"

$(5

55

$!

"*

$!E

!

!55!

其中&"!而&,!$L&JC(5J!"L("!M

$

1!J!"L"!%3;

点的形函数:如此可在有限的插值点上获得指数型的收敛速度

:

!%"!空间离散

首先!如同有限元方法!我们将整个求解区域分解成%&个互不重叠的四边形区域"&!"L+"&!我们将这种区域分解称之为四边形分解

(并用B标记!这里(代表区域分解后的网格参

数:相应的!测试函数被定义在单元"&上!并标

记为&("&

;为便于计算!需要引入参考单元:我们定义)L’P"!"(

!这样!参考单元可以用)%3来标记!当空间维数%3分别为$和!时!参考单元就分别为正方形和立方体:对于任意单元"&(B(!假定在单元局地坐标和参考单元间)

%3均存在一个可逆的转换函数G&Q)%3,"&:这样!我们就可以

通过这个变换完成实际物理区域和参考单元之间的转换!使计算统一在参考单元中进行:通过空间

分解B(!我们标记&(

!为定义在单元上的尝试函数空间!而’(

为定义在单元上的测试函数空间:让(!取和’(相同的函数基!就得到弱形式的离散化的P/1;AE03表示#

"&(!#’($N6"&(!!($L"&(!"$N"&(

!$$%!!">$其中!("!$(&(!!

&((’(

;区域分解以后!为提高近似解的精度!我们可以在单元区域上再增加节点"插值点$来增加近似解的阶次!这在有限元中被称为高阶元方法:如何选择增加节点的位置以及如何选择形函数将谱元法从普通的高阶有限元法中区分开来:为简便计!以下我们只考虑二维问题:

图"显示了基于M;4;3<A;展开的一个参考单元上的节点配置示意图!其中正交多项式的截断阶次为Dh+:可以看到!在采用M;4;3<A;展开的谱元法中!为了保持谱方法的收敛速度!并且增加数值计算的快捷!单元内节点"插值点$被放置在特定的配置点上!它们是由每一维上的M;4

;3<’;’P/9??’MCX/22C点通过张量积的形式构成的’$

$

(

:而每一维上的基函数由有限项正交多项式的和来

表示!再通过张量积的形式构成单元上任意配置

&A

!

"#

&"

%&

D

(

D

"&

5M(

"2!D*!+$

"&

5M5M

"2&"5M!*!+$#)*!+"&!!$DL#

(

D

(D

%&

D

&L"5!ML#

**&

.

$""""$"$

6&!!$D"DL

(

D

(D

5!&L"ML#

**

(

&D

"&

5M5M(

"%"QQ!D*!+$*!+$

这里)/=CX0/3矩&是每一单元进行坐标变换时的]阵!权重///55MLM!而

.

波传播问题!中心差分法被认为是最有效的直接积

,(

$分方法’;时间点"%N"!的解是通过解下面的线.

表示将原先本地单元上的

%3上实施#梯度算子在参考单元)

性方程组直接得到的!

.

"&

5M5M

"2)&"5M!*!+$*!+.0

R

&(D

"&

L&(

D

"&

".$

L

.

(.

&D

"".$2(G&

&.

$

!""!

$和有限元方法一样!通过单元聚合!我们最终得到了如""+

$式所表示的二阶线性常微分方程组:*+"!$N,

-"!$L."!$!""+

$这里!-"#$/3-#4!0-"#$/30-#4是相应的初值条件;而未知向量-包含了所有单元中所有M;4

;3<A;配置点上的离散解*-:*是全局质量矩阵!,是全局刚度矩阵!而.是力向量:来自那些单元共享节点的贡献在全局矩阵中被相加以满足单元边界上的连续性要求:在等参变换的帮助下!所有的单元矩阵计算可以统一在参考单元中实施:

在采用M;4;3<A;展开时!一个很重要的结果是由于使用了M;4;3<A;’P/9??’MCX/22C配置点技术并采用协调积分方式!得到的全局质量矩阵始终是个对角阵!这将导致我们可以直接得到线性方程组的解和采用显式的时间迭代方式!极大地节约了内存需求和计算时间:

%#!时间离散

为了求解""+$式所示的线性方程组!必须在’#!F(时间区间上进行积分:我们采用直接积分法!具体把整个时域离散成DF个时间点!!%L%.!!%%%DF!这里.!LF,DF;

由于采用M;4;3<A;展开时得到的全局质量矩阵*是对角阵!我们可以用中心差分法来实现""+$式的时间积分;中心差分法是一种条件稳定算法!然而!当*是对角阵时!它是显式的;当对时间步长的限制不是太严厉时!如同我们正在讨论的弹性

%N

".

!$

*$-LG%

P"

,P

%.!$$

-""&

$P%P

".

!$

$-;

%$!逐元技术

逐元概念在有限元方法中就有表述’,(

!D;A0/30在">>%年首先把它引入基于(G;XJ

?G;Z展开的谱元法中’"*(

:逐元技术主要是一种思想!在不同的方

法中可以有不同的应用形式!并产生不同的效果:本文首次将逐元技术引入基于M;4;3<A;展开的谱元法中!并作详细的阐述和效果分析:

传统的谱元法或有限元法解波动方程时!需要形成全局刚度矩阵和质量矩阵!考虑衰减时还需引进阻尼矩阵!不管是用一维等带宽存贮还是变带宽存贮!都需要耗费大量的存贮空间!尤其在谱元法中!随着阶次的提高!相应的带宽成倍增加:注意到在波动方程的时间迭代过程中!需要做多次矩阵和向量乘积"参看""&$式$!如122J!这里1和2J分别为全局刚度"或质量$矩阵和全局向量;考虑到全局矩阵的稀疏特性!全局的矩阵和相量相乘耗费了大量的计算时间;这些都极大地阻碍了谱元法在具体问题中的应用;我们注意到!在谱元法中!全局矩阵是单元矩阵按特定的结构聚合而成的;因此从单元的观点来看!122J的结果全局向量可以通过一系列单元层次的计算得到!

122JL

"*-1

&$2J

L

*"

-1&

22.&

J

$&

&

"",

$L

*

3/&

JL3J!

&

其中!-1&和2.&J分别为构成全局矩阵和全局向量的

单元矩阵和向量!而*代表聚合;在上述公式中!

&

!!#.

"#

&$

!

全局的矩阵向量乘积是在局部单元上完成的;具体步骤如下!首先!计算单元矩阵!而单元向量通过

&全局和局部的关系获得!也就是说!单元向量2.J&通过特定的关系从全局向量2.J扩散得到*然后!

也就是说!通过采用逐元技术!计算需求的减少和内存需求的减少同样巨大:

另外!由于在基于M;;3<A;展开的谱元法程序4中!绝大部分的工作可以局限在局部单元上实施!各单元间很容易获得负载平衡!并可通过聚合和扩散完成整个计算!因此采用逐元技术的谱元程序从结构上就具有优良的并行效率!具备大规模并行计算的特殊基础:

在局部单元上计算单元矩阵向量积!得到局部单元

&&&-22结果向量31/.JLJ*最后!局部的单元结果向量&3/J再通过特定的关系聚合成全局结果向量;就这样!

凭借在全局和局部之间正确的扩散和聚合!我们没有使用全局矩阵1就得到了需要的全局向量:因此!在程序中我们就没有必要集成并存贮全局矩阵!而只需计算并存贮单元矩阵:单元矩阵是稠密的!这极大地缩减了内存需求和提高了计算效率:

为说明这点!我们以一个二维正方形计算区域为例:在基于M;4

;3<A;展开的谱元法中!"被离散为7O7个互不重叠的四边形单元:假定M;4’3<A;展开的阶次为D!计算区域具有的总节点数

为"7DN"$O"7DN"$:这样!全局矩阵1的大小

为1$i"7DN"$O"7DN"$

($

<L’;而假如我们采用逐元技术!每个单元矩阵-1&的大小为’$i"

DN"$O"DN"$($

!7O7个单元的总需求为1&L7O7O’

$i"DN"$O"DN"$($

;两种方式的内存需求比如下式所示!

1&7$

O+DN"!""%

$&

$

假定7L"##!DL,!则1<

1&

&&+##:也就是说!在

此例中!通过采用逐元技术!内存需求至少降低了###倍:

现在让我们来关注计算需求!在基于M;4;3<A;展开的谱元法中!最大的计算负荷是计算全局的矩阵向量乘积122J!这里只关注耗时的乘法运算!在全局方法中需要进行的乘法总次数为F<L

’$i"7DN"$O"7DN"

$($

;采用逐元技术后!一个全局矩阵向量相乘被转换成7O7次单元矩阵向

量相乘!它们需要的总的乘法次数F<L7O7O

’$i"

DN"$O"DN"$($

;因此!两种方法的计算量之比大约为

F+<F&7$

!&ODN

"$

""*

$!数值算例

为了验证基于逐元技术的M;4;3<A;谱元法的精度和效率!我们选择了几个特定的计算模型!并将谱元法的结果和精确解或其他数值方法的结果进行对比:需要说明的是!由于采用了逐元技术!下述算例均在微型机上完成:

%!!均匀介质

首先!假定有一脉冲力源作用在弹性均匀介质中!对此情况有解析解:我们将之和谱元法计算的结果进行对比:计算模型见图$!它由一边长h$&,#.的正方形区域构成!以区域中心为坐标原点!边界条件为放射边界!其纵)横波速度和密度分别如下!4.h

$>##.,?!4+h","".,?!h"

>##E4,.!

:一个具有有限带宽的水平力作用在坐标为"b!##.!b!##.$的点#上:力源的波

形函数为P/9??函数的一阶导数!其中心频率为

#‘7!我们将记录点8的坐标定为"!##.!!##.$

:在谱元法中!计算模型被离散成,+i,+个正方形单元!每个单元的尺寸为+#.i+#.:在每个单元上!采用M;4;3<A;正交多项式的阶次为L+:而时间步长为.!h#:*.?:吸收带被加在计算区域的每一个边界上用来模拟放射性边界条件:

图!显示了采用谱元法计算时!水平速度分量SC在#:!,?时刻的声场快照图!在那里我们可以清晰地辨认出纵波和横波:图+"/$和图+"X$分别比较了应用解析算法和谱元法两种方法在接收点8得到的水平速度分量和垂直速度分量的波形:在每个图上!实线标记解析解的结果!虚线标记谱元法得到的结果!而点线标记两者之间的误差!

";""B#&"D

!

"#

&!

并且为了显示方便!误差被放大了!倍:从图+中可以看到!谱元法的结果和解析方法得到的精确

解符合得很好!对于整个记录波形!最大相对误差不超过$

[:

图"!

均匀介质模型

图#!在,%#&>时刻45

快照图

图$!对于水平力源在弹性均匀介质中激发的声场&谱元法解和解析解比较

"$水平速度分量*"$垂直速度分量/X

"%"!带有一水平分层的分区均匀介质

现在我们考虑带有一水平分层的分区均匀介质情况:计算模型显示在图&中!它仍然为一边长为

8$的坐标为":&##.!$%&.$

这里!伪谱法"被用来和谱元法计算的结果$OD作比较:在计算空间微分时!伪谱法在水平和垂直轴上均采用NC9A0;A变换:对于图&所示的规范的分区均匀介质模型!NC9A0;A变换被认为在直至每最小波长$个采样点$下均具有I90?2采样率"Jf

无限精度"机器精度$:因此可以认为伪谱法在这个模型下的结果是精确的:

在伪谱法中!整个计算区域被划分为$&,i$&,个正方形网格!每个网格的边长为"#.:时间步长

Bh$&,#.的正方形区域!以区域中心为坐标原点!

边界条件为放射边界!但在区域中存在一水平分层将计算区域一分为二:下方地层为慢速地层!其纵)横波速度及密度分别为4,!,!$>##.","".?4?"h+"h.

!,">##E.:上方地层为快速地层!其纵)横波速4"h#

度及密度分别为4,!4,!,##.$#&%.??$h!+$h.

!,$,*#E.:算例$:"中的水平带限声源仍然被加载4$h#

到位于慢速地层中#点上!其坐标为":#.!b$$&.$我们使用分别位于慢)快地层中的两个记录点来记录接收波形:记录点8!记录点&##.!b$%&.$"的坐标为"

!h#:*.?:为了避免NC9A0;A变换带来的周期环.

叠!在水平和垂直边界上都加上了吸收带:

在谱元法中!整个计算区域被离散成,+i,+

"#

&+

!

个正方形单元!在每个单元上!采用M;;3<A;正4交多项式的阶次为DL+:考虑到在伪谱法和谱元法中对分层边界的位置定义不同!伪谱法定义在网格中间!而谱元法则定义在单元边界上:在这里!谱元每个单元的尺寸被规定为&#.i&#.:这样!谱

元法的计算区域将大于伪谱的计算区域!但由于声源和接收器的相对位置不变并采用了放射边界!这不会影响我们的比较:谱元法中的时间步长也是

!h#:*.?:同样!吸收带被加在边界上用来模拟.

放射性边界条件

:

图0!

带有一水平分层的分区均匀介质模型

图&!在,%#&>时刻45快照图%’72-’

:!,?时刻水平!!图,显示了由谱元法计算的在#

速度分量SC的声场快照图!在图中我们可以清楚地分辨出直达纵波)横波!反射纵波和横波!折射纵波"$"$和横波以及各种头波:图%和图%分别比较了伪/X

谱法和谱元法两种方法在接收点8"得到的水平速度分量和垂直速度分量的波形:图*和图*则分别"$"$/X

比较了伪谱法和谱元法两种方法在接收点8$得到的水平速度分量和垂直速度分量的波形:在每个图上!实线标记伪谱法解的结果!虚线标记谱元法得到的结果:所有比较都获得了非常一致的结果:可见谱元法仅使用较低阶的正交多项式展开就在这种分区均匀介质上获得了和伪谱法相当的精度

:

图?!对于水平力源在带水平分层的分区均匀介质中激发&

6!点接收波形&谱元法解和伪谱解比较

水平速度分量*"垂直速度分量"$$/X

!

"#

&&

图@!对于水平力源在带水平分层的分区均匀介质中激发&

6"点接收波形&谱元法解和伪谱解比较

"$水平速度分量*"$垂直速度分量/X

"%#!带有一倾斜分层的分区均匀介质

这里我们考察带有一倾斜分层的分区弹性均匀介质情况:计算模型显示在图>中!它仍然为一边长为Bh$&,#.的正方形区域!以区域中心为坐标原点:但在区域中存在一倾斜分层将计算区域一分为二!分层通过区域中心!和水平轴形成的倾角为"*:+!d:下方地层为慢速地层!上方地层为快速地层!它们各自的纵)横波速度及密度均和例$:$相同:算例$:"中的水平带限声源仍然被加载到位于地层中#点上!其坐标为"#.!b!##.$

:

样!它在谱元法中是自然满足的:

图"#分别显示了计算得到的在时刻

#:"&!#:$&!#:!&和#:+&?的水平速度分量SC的声场快照图:弹性波在这种带有一倾斜分层的分区均匀介质中的传播过程被清楚演示!这也显示谱元法对于模拟复杂结构介质中的弹性波传播具有独特优势

:

图A!带有一倾斜分层的分区均匀介质模型

由于存在倾斜分层!只有谱元法被用于计算此模型:这里!整个计算区域被离散成,+i,+个任意四边形单元!而不再是正方形单元:在每个单元区域上!我们仍然采用阶次为DL+的M;;3<A;正交4多项式来进行插值:时间步长取.!h#:&.?:同时

自由边界条件被用于本计算模型中!和有限元法一

"%$!具有弯曲表面的均匀介质

最后!我们将谱元法用于具有起伏地表的实际

图!,!45声场快照图

"$$$$*"*"*"/X=<#:"&?#:$&?#:!&?#:+&?

"#

&,

!

地层中的弹性波传播模拟中!这里只考虑二维情况:图""所示为具有一平缓小山的二维地形剖面图!模型的尺寸为+###.i",##.!而小山的高度为",#.!整个地表由一个P/9??函数描述:地层均匀!其纵)横波速度和密度分别如下

!

,!,!,’?’?.:计!$##."*+%:&.$##E4+h.L#h$

算模型被离散成&#i$#个单元!其网格划分见图

#!结论

谱元法通过在每一个单元上应用谱展开从而结合了谱方法和有限元方法的优点:在本文中!详细阐述了基于M;;3<A;正交多项式展开的谱元法基本4理论及相应数学公式:由于采用了协调积分和M’;4;3<A;’P/9??’MCXC22C配置点技术!得到了对角形式

的全局质量阵!简化了运算:逐元技术被首次用于"":在每个单元上!采用M;;3<A;正交多项式的阶4

次为DL,:而时间步长为.!h#:&.?:吸收带被加在除上表面外的每一个边界上用来模拟放射性边界条件:算例$:"中的水平带限声源被加载到位于地层中#点上!中心频率为%‘7其坐标为""!,#.!++#.$

:

图!!!具有起伏地表的二维地形剖面及其网格划分图

图"$分别显示了计算得到的在时刻#:"!#:$!

:!和#:+?的水平速度分量SC的声场快照图:弹性波在这种具有起伏地表的地层中的传播过程被清楚演示!这也显示谱元法完全胜任在具有弯曲界面的复杂介质中的传播模拟工作

:

图!"!45声场快照图

"/$#:"?*"X$#:$?*"=$#:!?*"<$#:+?

基于M;4;3<A;多项式展开的谱元法程序中!并进行了详细阐述和说明!它极大地减少了内存和计算需求!对谱元法的推广应用具有重大意义:最后!通过四个数值算例验证了我们这种谱元法的高精度和强适应性!并由此证明这种应用逐元技术的谱元法适合用于大型复杂介质和结构中的弹性波传播模拟:

致谢!本文的完成得到了意大利国立海洋和地球物理研究所P;7/D;A0/30博士的无私帮助!并就谱元法给予很多建议和讨论!在此作者表示真诚感谢:

附录

M;4;3<A;正交多项式3BJ"

C$!JL#!"!14是如下奇异性D29A.’M0C9Z011;方程的特征函数’$

$

(

"""PC$

$BTJ"C$$TNJ"JN"$BJ"C$L#!")"

$注意到它的权函数U"C$L"!这使它在进行数值积分时很有利:如将BJ""$进行归一使得BJ""$L"!这样J阶M;4;3<A;多项式可以如下定义

’J,$

(

BJ"

C$L$J*"P"$"JCP

!")$

$"L#"$"

"

$JP$"J$"!J

$

其中!’J,$(表示对J,$取整:而M;4;3<A;’P/9??’MCX/22C点及相应的权重定义如下

C#LP"

!!CDL"!C1!D"P"$取BT)!

$M"

ML"!D的零点!/ML

D"DN"$’B!!"1!D$!")+$D"CM$

($

ML#!其中BTD分别为D阶M;4;3<A;多项式的一阶导数:而一维上的基函数,5".$

可以如下定义!"#

!

,5

".$L"""DD$

$BTD"$N"B.;")&$D5.P.

5参!考!文!献

"!‘0A?=G(:I9.;A0=/1(C.W

92/20C3CHS32;A3/1/3<-52;A3/1N1C6?:YC1:"!I;6aCAE#]CG3B01;J^?

C3?!">**$!Y0A0;95]:O’DY6/Z;WACW/4/20C303G;2;AC4

;3;C9?.;<0/#Y;’1C=02J’?2A;??H0302;’<0HH;A;3=;.;2GC<:P;CWGJ?0=?!">*,!&"#**>+>#"

!!M;Z/3<;A)Q:NC9A2G’CA<;AH0302;’<0HH;A;3=;O’DY?;0?.C’

4A/.?:P;CWGJ

?0=?!">**!&!#"+$&+"+!,+!王秀明!张海澜:用于具有不规则起伏自由表面的介质中弹性

波模拟的有限差分算法:中国科学!P辑!$##+!!+"&$#+*"++>!

&!N031/J?C3L)!D=A0Z;3M-:TG;.;2GC<CH6;04

G2;<A;?0<9/1?’/A;Z0;6:)WW

1R;=G!">,,!Q;Z">#%!&+%+*,!‘94

G;?T]Q:TG;N0302;-1;.;32R;2GC<!M03;/AD2/20=/3<@J3/.0=N0302;-1;.;32)3/1J?0?:-341;6CC<(10HH?!I]#OA;3’20=;’‘/11S32;A3/20C3/1!">*%

%!R9aP:-1/?20=6/Z;.04

A/20C3602GH0302;;1;.;32.;2GC<:)=2/P;CWGJ

?0=/D030=/!">*+!$%#$,*+$%**!T;34a(:TGA;;’<0.;3?0C3/1H0302;;1;.;32/3/1J

?0?CH6/Z;?03/3/=C9?20=.;<0/602G03=19?0C3:])=C9?2DC=).!">**!*,#+"+++$$

>!PA/Z;?Q:D0.91/20C3?;0?.0=6/Z;WACW/4

/20C303!@;1/?20=.;<0/9?034?2/44;A;<’4A0<H0302;<0HH;A;3=;?:L911D;0?.DC=).!">>,!*,#"#>"+""#,

#!R/AH9A2F]:)==9A/=JC

HH0302;’<0HH;A;3=;/3<H0302;’;1;.;32.C<;11034CH2G;?=/1/A/3<;1/?20=6/Z;;f9/20C3?:P;CWGJ?0=?!">*+!+>#&!!+&+>

"!FC./202?=G@!Y01C22;]O:TG;?W

;=2A/1;1;.;32.;2GC<#)3;H’"#

&%

H0=0;322CC12C?0.91/2;2G;?;0?.0=A;?WC3?;CH$@/3<!@4;C1C4’0=/1?2A9=29A;?:L911D;0?.DC=).!">>*!**#!,*+!>$"$!UA?7/4D):DW;=2A/1.;2GC<?HCAWACX1;.?03=C.W

1;54;C.;’2A0;?:](C.W92OGJ

?!">*#!!%#%#+>$"!!FC?1CHH@!F;??1;A@!N01GC)g!;2/1:DC1920C3CH2G;;f

9/’20C3?CH<J3/.0=;1/?20=02JXJ/(G;XJ=G;Z?W;=2A/1.;2GC<:P;C’WGJ

?0=?!">>#!&&#%!++%+*"+!(/A=0C3;]R!B/34O]

:)(G;X?G;Z=C11C=/20C3.;2GC<HCA2G;6/Z;;f9/20C3034;3;A/107;<=CCA<03/2;?:(C.W92;A?/3<N190<@J

3/.0=?]C9A3/1!">>!!$#$,>+$>#"&!(/A=0C3;]R:@C./03<;=C.WC?020C3HCA6/Z;WACW/4

/20C3WACX’1;.?:]C9A3/1CHD=0;320H0=(C.W92034!">>"!,"+$#+&!++%$",!O/2;A/)T:)?W;=2A/1;1;.;32.;2GC<HCAH190<<J

3/.0=?#M/.03/AH1C603/=G/33;1;5W/3?0C3:]CH(C.W92OGJ?0=?!">*+!&+#+,*++**

"%!D;A0/30P!OA0C1C-:)39.;A0=/103Z;?204/20C3CH(G;XJ

?G;Z?W;=2A/1;1;.;32.;2GC<HCA/=C9?20=6/Z;WACW/4

/20C3:OAC=CH"!2G

SR)(DBCA1<(C34A;??C3(C.W

92/20C3/3<)WW10;<R/2G;’./20=?:@9X103!SA;1/3<:">>"

"*!D;A0/30P:)W/A/11;1?W

;=2A/1;1;.;32.;2GC<HCA/=C9?20=6/Z;.C<;1034:]CH(C.W

2)=C9?20=?!">>%!&""$#&!+,>">!D;A0/30P:@C9X1;’4A0<(G;X?G;Z?W

;=2A/1;1;.;32?HCA/=C9?20=6/Z;.C<;1034

:B/Z;RC20C3!$##+!!>#!&"+!,#$#!@/9E?G;AB!-.;AJ)N:)==9A/=J03.C<;110342

G;/=C9?20=6/Z;;f9/20C3602G(G;XJ?G;Z?W;=2A/1H0302;;1;.;32?:N0302;-1;.;32?03)3/1J?0?/3<@;?04

3!">>%!$,#""&+"$*$"!FC./202?=G@!TAC.W]:S32AC<9=20C32C2G;?W

;=2A/1;1;.;32.;2GC<HCA2GA;;’<0.;3?0C3/1?;0?.0=6/Z;WACW/4/20C3:P;CWGJ?]S32!">>>!"!>#*#,+*$$

$$!(/392C(:DW;=2A/1R;2GC<?03N190<@J

3/.0=?:I;6aCAE#DWA034;A’Y;A1/4

!">*,""

本文来源：https://www.bwwdw.com/article/znhq.html