求解弹性力学方程的主要方法及其特点比较

求解弹性力学有类方程，共15个方程。3个平衡方程，6个物理方程，6个几何方程。

弹性力学是固体力学的重要分支，它研究弹性物体在外力和其他外界因素作用下产生的变形和内力，又称弹性理论。它是材料力学、结构力学、塑性力学和某些交叉学科的基础，广泛应用于建筑、机械、化工、航天等工程领域。

弹性体是变形体的一种，它的特征为：在外力作用下物体变形，当外力不超过某一限度时，除去外力后物体即恢复原状。绝对弹性体是不存在的。物体在外力除去后的残余变形很小时，一般就把它当作弹性体处理。

弹性力学所依据的基本规律有三个：变形连续规律、应力-应变关系和运动(或平衡)规律，它们有时被称为弹性力学三大基本规律。弹性力学中许多定理、公式和结论等，都可以从三大基本规律推导出来。

一、变形连续规律 弹性力学（和刚体的力学理论不同）考虑到物体的变形，但只限于考虑原来连续、变形后仍为连续的物体，在变形过程中，物体不产生新的不连续面。如果物体中本来就有裂纹，则弹性力学只考虑裂纹不扩展的情况。

反映变形连续规律的数学方程有两类：几何方程和位移边界条件。几何方程反映应变和位移的联系，它的力学含义是，应变完全由连续的位移所引起，在笛卡儿坐标系中，几何方程为：

若所考虑的物体Q在其一部分边界B1上和另一物体Q1相连接，而且Q在B1上的位移为已知量，在B1上便有位移边界条件:

二、应力-应变关系 弹性体中一点的应力状态和应变状态之间存在着一定的联系，这种联系与如何达到这种应力状态和应变状态的过程无关，即应力和应变之间存在一一对应的关系。若应力和应变呈线性关系，这个关系便叫作广义胡克定律，各向同性材料的广义胡克定律有两种常用的数学形式：

和

式中为应力分量；λ和G为拉梅常数，G又称剪切模量；E为杨氏模量（或弹性模量）；v为泊松比（见材料的力学性能）。λ、G、E和v四个常数之间存在下列联系：

三、运动（或平衡）规律 处于运动（或平衡）状态的物体，其中任一部分都遵守力学中的运动（或平衡）规律，即牛顿运动三定律，反映这个规律的数学方程有两类：运动（或平衡）微分方程和载荷边界条件。在笛卡儿坐标系中，运动（或平衡）微分方程为：

对于均匀而且各向同性的物体，应力分量可按式（3a）用应变分量表示，而应变分量又可按式（1）用位移分量表示。两个公式依次代入方程（5），便得到用位移表示的运动微分方程：

式中θ为体应变，即:

△为拉普拉斯算符，即:

类似地，在方程（6）中略去惯性力，便可得到用位移分量表示的平衡微分方程。

如果考虑物体一部分边界B2是自由的，在它的上面有给定的外载荷，则根据作用力和反作用力大小相等方向相反的原理，在B2上有如下载荷边界条件：

对弹性力学的动力问题，还需说明物体的初始状态，即:

当t=t0时，

柱体扭转和弯曲 一个侧面不受外力的细长柱体，在两端面上的外力作用下会产生扭转和弯曲。根据圣维南原理，柱体中间部分的应力状态只与作用在端面上载荷的合力和合力矩有关，而与载荷的具体分布无关。

杨顶辉

（石油大学地球科学系，北京　100083）

滕吉文　张中杰

（中国科学院地球物理所，北京　100101）

摘要　快速有效地模拟地震波在各向异性介质中的传播在现今勘探地震学中具有重要的意义。一种算法的稳定性分析对于加快计算速度非常必要。本文首先利用矩阵和向量来描述波传播方程，针对二维和三维一般各向异性介质中的弹性波方程，提出了一种快速且占用内存少的有限差分方法；然后系统地研究了二维均匀、非均匀各向异性情况下波动方程有限差分格式的稳定性条件；进一步给出了某些特殊各向异性情况下有限差分方法的稳定性具体公式。最后，本文也对三维有限差分格式的稳定性问题进行了研究。

关键词　弹性波方程　有限差分　稳定性　各向异性介质

1　引言

地震波传播的数值模拟在地球科学中具有重要的意义。在各向异性地震模拟的各种方法中，基于Kennett研究工作[11,12]的反射率方法是最流行的数值技术之一。基于走时方程渐近解的射线追踪方法是模拟地震各向异性的另一种有效方法[5,6]。Kosloff等[13]利用Fourier方法模拟了地震各向异性。而Chen[7]则使用有限元方法模拟非均匀各向异性问题。虽然有限差分方法已被广泛应用于各向同性介质中的弹性波模拟，但是利用这种方法来模拟地震各向异性问题并不普遍。Tsingas等[16]利用有限差分算子发展了一种模拟算法以求解横向各向同性介质中的偏微分方程，这种算法是基于MacCormack型的分离格式[2]。Faria等[8]基于交错网格格式[17]，利用有限差分算法模拟了二维横向各向同性问题。最近，Igel等[9]基于褶积算法给出了一种模拟一般地震各向异性的有限差分算法。

快速和少存贮量是有限差分方法的优点。随着大尺度波场模拟的需要和大规模并行计算的发展，复杂介质或高维（二维和三维）模型的有限差分地震模拟不可避免。虚谱法因其空间算子准确地达到Nyquist频率而深受欢迎，然而虚谱法需要富利叶变换，从而对于三维各向异性模拟非常耗时。同时，采用富利叶变换意味着每一个网格点与其它的点相互影响。在某种意义上，这与动力学的局部弹性性质不一致。因此当我们为地震模拟设计有限差分格式时，考虑差分算子的局部性是必要的。另一方面，考虑差分算子的局部性对于提高算法的并行性非常重要。因为最邻近网格点间的信息交换最快，因而对于各向异性大尺度模型波场的模拟是可行的。

基于上述原因，本文针对一般各向异性介质中的地震波传播问题，给出了一种快速且占用内存少的有限差分方法。事实上，这是算法[10,18]的一种推广。

通常地，时间推进算法被使用于地震波传播的数值模拟中，为了保证算法稳定，时间增量受算法稳定性条件的限制。选择合适的时间步长不仅能保证数值计算稳定，而且能加快计算速度。否则的话，不但会产生非物理数值振荡，甚至会导致错误的结果。合理时间增量的选取决定于差分格式和描述介质特征的介质参数或速度，换句话说，决定于差分格式的稳定性条件。因此有限差分格式的稳定性问题在数值计算中十分重要。尽管这一问题在文献[10,18]中针对某些特殊情况作过研究，但他们并没有对一般各向异性问题的有限差分格式及其稳定性做过详细、系统的研究。我们的目的就是针对这一问题给出一般的有限差分格式及其稳定性条件，进一步地给出某些特殊各向异性情况下的稳定性公式。结果表明在各向同性情况下我们的结果与Aboudi[1]的结果是一致的。

2　有限差分格式

21　三维各向异性

各向异性介质中的运动平衡方程可成如下形式：

岩石圈构造和深部作用

其中：ρ是介质密度；fx,fy和fz分别表示力源在x,y和z方向上的分量；ux,uy和uz分别表示x,y和z方向上的位移分量。

应力应变关系为 ，其中弹性参数矩阵 ，并且ci,j=cj,i,i,j=1,2，…，6；

应力矩阵σ＝（σxx，σyy，σzz，σyz，σxz，σxy）T；

应变矩阵ε＝（εxx，εyy，εzz，εyz，εxz，εxy）T；

且应力与位移的关系为：

岩石圈构造和深部作用

令

岩石圈构造和深部作用

岩石圈构造和深部作用

那么方程（1）可写为：

岩石圈构造和深部作用

显然，A,E,Q,B+D,C+G和F+H是实的对称矩阵。在不考虑源项 的前提下，采用有限差分逼近方程（2），可得下列有限差分格式：

岩石圈构造和深部作用

其中：Δx，△y和△z分别表示x,y和z方向的空间步长，△t表示时间步长，并且，

岩石圈构造和深部作用

22　二维各向异性

类似地，二维各向异性介质中的波动方程可写为：

岩石圈构造和深部作用

并且有下列差分格式：

岩石圈构造和深部作用

3　稳定性条件

31　均匀介质中二维有限差分格式的稳定性条件

均匀介质情况下格式（5）可简化为：

岩石圈构造和深部作用

根据Richtmyer和Morton的稳定性理论分析，我们令

岩石圈构造和深部作用

其中U0是一常数向量，方程（6）变为：]]

2I—2Pλ＋I）U0=0

其中I表示一单位矩阵，

岩石圈构造和深部作用

其中： ， ， ， ，

若系数行列式为0，则满足：

岩石圈构造和深部作用

定理1　差分格式（6）稳定的条件是：

岩石圈构造和深部作用

其中函数F（α，β）和α，β为：

岩石圈构造和深部作用

岩石圈构造和深部作用

其中k=Δz/Δx

证明：据差分格式（6）的稳定性，方程（7）中的λ满足‖λ‖≤1。

根据（7）和引理2（见附录），有下列不等式：

岩石圈构造和深部作用

由A、Q和C+G的对称性，可知矩阵 也为对称的，则由引理3（见附录）和不等式（8），有

岩石圈构造和深部作用

由矩阵A，Q和C+G的元素的非负性，可令0≤α， ，则要使（9）成立，只须

岩石圈构造和深部作用

令

岩石圈构造和深部作用

并令偏导数

岩石圈构造和深部作用

根据波动方程的特性，有

‖A‖＞0,‖Q‖＞0,‖C+G‖＞0

所以（0，0）和（ ， ）可能为函数的极值点。显然，

f(0,0)=0

岩石圈构造和深部作用

下面我们讨论（α，β）≠（0，0）或（ ， ）的情况：

由（11）和（12）式，如果 ，那么

岩石圈构造和深部作用

其中α，β可能是f（α，β）的极值点，故稳定性条件为：

岩石圈构造和深部作用

否则 是函数的最大值点，且稳定性条件为：

岩石圈构造和深部作用

特别地，当Δx﹦Δz时，有下列简化的稳定性条件：

岩石圈构造和深部作用

其中α和β由（13）和（14）决定，且k=1。

32　非均匀情况下的稳定性条件

在非均匀情况下，对于差分格式（5）的稳定性条件是难以确定的。然而根据偏微分方程的数值方法理论[15]，我们可以采用“冻结系数”法[14]来分析其稳定性。进一步地，我们给出非均匀介质情况下差分格式（5）的稳定性条件。

事实上，如果介质参数函数为连续有界的，则通常我们可以近似地将一个小的计算区域看成是均匀的，那么差分格式（5）可以退化成格式（6），这样在小区域范围内，我们可以像均匀情况一样获得稳定性条件。进一步根据介质参数函数的连续有界特性，我们可以获得差分格式（5）的稳定性条件为：

如果 那么

岩石圈构造和深部作用

否则，

岩石圈构造和深部作用

其中H（α，β）、α和β为：

岩石圈构造和深部作用

33　某些特殊介质中差分格式的稳定性条件

显然，通过计算格式（6）中矩阵A、Q和C+G的范数，可以分别地获得各向同性和横向各向同性情况下的稳定性条件为：

各向同性

岩石圈构造和深部作用

或

岩石圈构造和深部作用

其中λ，μ为拉梅常数， 是P波速度。

显然，稳定性条件（15）与Aboudi[1]的结果一致。

横向各向同性

如果 那么 ≤p，否则

岩石圈构造和深部作用

其中，

岩石圈构造和深部作用

其中A,N,L,F和C为弹性常数。

类似地，我们可以获得非均匀和其它特殊各向异性介质（如：立方体各向异性、正交各向异性等介质）情况下的稳定性条件。

4　三维各向异性情况下差分格式的稳定性条件

像前面二维情况一样分析可得如下稳定性条件：

定理2　如果

max[k1·‖A‖＋k2·‖E‖＋k2·‖Q‖，f2（θ2，θ2，θ3）]≤1，那么均匀介质情况下差分格式（3）是稳定的。其中 ，函数f2（θ1，θ2，θ3）被定义为：

岩石圈构造和深部作用

其中： ， ， ， ， · ， ；

并且θ1，θ2，θ3满足：

岩石圈构造和深部作用

其中：a=4k1·‖A‖，b=4k2·‖E‖，d=4k3·‖Q‖，g=4k4·‖B+D‖，e＝4k5·‖C+G‖，f=4k6·‖F+H‖。

对于三维非均匀介质情况，经由“冻结系数”法可类似地分析其稳定性条件。

5　总结和讨论

数值模拟地震波传播的有限差分方法是一种重要的工具，而其稳定性条件是提高计算速度的关键之一。然而对于一般二维和三维各向异性介质情况，系统深入地研究其差分格式和稳定性条件尚少，本文给出了一种快速且占有内存少的有限差分格式，并系统地分析和推导了一般均匀和非均匀各向异性情况下差分格式的稳定性条件。我们相信本文获得的结果有助于各向异性数值模拟的发展，并为有限差分方法的广泛应用提供理论依据。

参考文献

[1]JAboudiNumerical simulation of seismic sourcesGeophysics,1971,36,810～821

[2]ABayliss,KEJordan,BJLeMesurier and ETurkelA fourth-order accurate finite-difference scheme for the computation of elastic wavesBullSeisSocAm,1986,76,1115～1132

[3]DCBooth and SCrampinThe anisotropic reflectivity technique:theoryGeophysJRAstrSoc,1983（a）,72,755～756

[4]DCBooth and SCrampinThe anisotropic reflectivity technique:anomalous arrivels from an anisotropic upper mantleGeophysJRAstrSoc,1983（b）,72,767～782

[5]VCervenySeismic rays and ray intensities in inhomogeneous anisotropic mediaGeophysJRoyAstrSoc,1972,291～13

[6]VCerveny and PFirbasNumerical modelling and inversion of travel-times of seismic body waves in inhomogeneous anisotropic mediaGeophysJRoyAstrSoc,1984,76,41～51

[7]KHChenPropagating numerical model of elastic wave in anisotropic inhomogeneous media:finite element methodsThe 54th SEG Annual Meeting Expanded Abstracts,Houston,USA,1984,631～632

[8]ELFaria and PLStoffaFinite-difference modelling in transversely isotropic mediaGeophysics,1994,59,282～289

[9]HIgel,PMora and BRiolletAnisotropic wave propagation through finite-difference gridsGeophysics,1995,60,1203～1216

[10]KRKelly,RWWard,STreitel and RMAlfordSynthetic seismograms:A finite-difference approachGeophysics,1976,41,2～27

[11]BLNKennettTheoretical reflection seismograms for an elastic mediumGeophysProsp,1979,27,301～321

[12]BLNKennettSeismic wave propagation in stratified mediaCambridge:Cambridge University Press,1983

[13]DKosloff and EBaysalForward modelling by a Fourier methodGeophysics,1989,47,1402～1412

[14]陆金甫，顾丽珍，陈景良偏微分方程差分方法北京：高等教育出版社，1988

[15]RDRichtmyer and KWMortonDifference methods for initial value problemsNew York:Interscience,1967

[16]CTsingas,AVafidis and ERKanasewichElastic wave propagation in transversely isotropic media using finite differencesGeophysProsp,1990,38,933～949

[17]JVirieusP-SV wave propagation in heterogeneous media:Velocity-stress finite-difference methodGeophysics,1986,51,889～901

[18]ZJZhang,QDHe,JWTeng,et alSimulation of 3-component seismic records in a 2-dimensional transversely isotropic media with finite-differenceCanJExplGeophys,1993,29,51～56

附录：引理

引理1　对于实系数方程λ2—2dλ＋1=0,｜d｜≤1是它的根λ满足∣λ｜≤1的充分必要条件。

引理1的证明参见文献[14]。

引理2　若A∈Rn×n且A=A′，I是一个单位矩阵，则对于下列方程

∣λ2I—2Aλ＋I∣＝0，

‖A‖≤1是方程的根λ满足∣λ｜≤1的充分必要条件。

证明：充分性

因为A为实对称矩阵，存在T-1、T∈Rn×n使得

T-1AT=diag（d1,d2，…，dn）

根据引理2中的方程及上述等式，可获得

（λ2—2d1λ＋1）…（λ2—2dnλ＋1）=0

由‖A‖≤1和ρ（A）≤‖A‖，有‖ρ（A）‖≤1

根据 ，有

故对满足方程的任一根λ有∣λ∣≤1。

必要条件：因∣λ∣≤1且A为实对称矩阵，故由引理1可获得：

∣di∣≤1,i=1,2，…，n

即ρ（A）≤1。

又因A为正矩阵，所以ρ（A）＝‖A‖，即‖A‖≤1。

引理3　如果A∈Rn×n，并且A=A′，那么

（i）如果‖I—A‖≤1，则‖A‖≤2;

（ii）如果‖A‖≥0，则‖A‖≤2是‖I—A‖≤1成立的充分必要条件。

证明：（i）因为A为实对称矩阵，故存在T-1、T∈Rn×n使得T-1AT=diag（d1,d2，…，dn）≡D

根据‖I—A‖≤1，有‖T（I—D）T-1‖≤1

显然，I—D为一正规矩阵，所以‖T（I—D）T-1‖＝‖I-D‖≤1

即

所以， ，即‖D‖≤2。

又因为A为正规矩阵，所以

‖A‖＝‖TDT-1‖＝‖D‖，即‖A‖≤2

（ii）必要条件已在（i）中证明，下面证明充分条件。

由（i）的证明过程可知：

‖A‖＝‖D‖

因为‖A‖≤2并且‖A‖≥0，有0≤‖D‖≤2，即max∣di｜≤2

所以我们可获得 。

即：‖I—D‖＝‖T-1（I—A）T‖≤1

由A的对称性可得‖I—A‖≤1。

求解弹性力学有类方程，共15个方程。3个平衡方程，6个物理方程，6个几何方程。

弹性力学是固体力学的重要分支，它研究弹性物体在外力和其他外界因素作用下产生的变形和内力，又称弹性理论。它是材料力学、结构力学、塑性力学和某些交叉学科的基础，广泛应用于建筑、机械、化工、航天等工程领域。

弹性体是变形体的一种，它的特征为：在外力作用下物体变形，当外力不超过某一限度时，除去外力后物体即恢复原状。绝对弹性体是不存在的。物体在外力除去后的残余变形很小时，一般就把它当作弹性体处理。

弹性力学所依据的基本规律有三个：变形连续规律、应力-应变关系和运动(或平衡)规律，它们有时被称为弹性力学三大基本规律。弹性力学中许多定理、公式和结论等，都可以从三大基本规律推导出来。

一、变形连续规律 弹性力学（和刚体的力学理论不同）考虑到物体的变形，但只限于考虑原来连续、变形后仍为连续的物体，在变形过程中，物体不产生新的不连续面。如果物体中本来就有裂纹，则弹性力学只考虑裂纹不扩展的情况。

反映变形连续规律的数学方程有两类：几何方程和位移边界条件。几何方程反映应变和位移的联系，它的力学含义是，应变完全由连续的位移所引起，在笛卡儿坐标系中，几何方程为：

若所考虑的物体Q在其一部分边界B1上和另一物体Q1相连接，而且Q在B1上的位移为已知量，在B1上便有位移边界条件:

二、应力-应变关系 弹性体中一点的应力状态和应变状态之间存在着一定的联系，这种联系与如何达到这种应力状态和应变状态的过程无关，即应力和应变之间存在一一对应的关系。若应力和应变呈线性关系，这个关系便叫作广义胡克定律，各向同性材料的广义胡克定律有两种常用的数学形式：

和

式中为应力分量；λ和G为拉梅常数，G又称剪切模量；E为杨氏模量（或弹性模量）；v为泊松比（见材料的力学性能）。λ、G、E和v四个常数之间存在下列联系：

三、运动（或平衡）规律 处于运动（或平衡）状态的物体，其中任一部分都遵守力学中的运动（或平衡）规律，即牛顿运动三定律，反映这个规律的数学方程有两类：运动（或平衡）微分方程和载荷边界条件。在笛卡儿坐标系中，运动（或平衡）微分方程为：

对于均匀而且各向同性的物体，应力分量可按式（3a）用应变分量表示，而应变分量又可按式（1）用位移分量表示。两个公式依次代入方程（5），便得到用位移表示的运动微分方程：

式中θ为体应变，即:

△为拉普拉斯算符，即:

类似地，在方程（6）中略去惯性力，便可得到用位移分量表示的平衡微分方程。

如果考虑物体一部分边界B2是自由的，在它的上面有给定的外载荷，则根据作用力和反作用力大小相等方向相反的原理，在B2上有如下载荷边界条件：

对弹性力学的动力问题，还需说明物体的初始状态，即:

当t=t0时，

柱体扭转和弯曲 一个侧面不受外力的细长柱体，在两端面上的外力作用下会产生扭转和弯曲。根据圣维南原理，柱体中间部分的应力状态只与作用在端面上载荷的合力和合力矩有关，而与载荷的具体分布无关。因此，柱体中间部分的应力有以下的表达式：