王圣柱 周建科 張奎華 王藝豪

1） 中國山東青島266580中國石油大學（華東）地球科學與技術學院2） 中國山東東營257061中石化勝利油田西部新區(qū)研究院 3） 中國武漢430079武漢大學測繪學院

?

有限元三角網(wǎng)格中波動的頻散分析及數(shù)值仿真

1） 中國山東青島266580中國石油大學（華東）地球科學與技術學院2） 中國山東東營257061中石化勝利油田西部新區(qū)研究院 3） 中國武漢430079武漢大學測繪學院

波動問題有限元離散后會引起數(shù)值誤差， 數(shù)值頻散的本質就是數(shù)值誤差傳播引起的非物理解． 數(shù)值頻散不僅沒有實際意義， 而且還會影響對真實波動現(xiàn)象的認識． 為厘清有限元三角網(wǎng)格中波動數(shù)值頻散的影響因素， 本文推導了集中質量矩陣和一致質量矩陣的頻散函數(shù)， 同時給出了組合質量矩陣的頻散函數(shù)， 并對不同質量矩陣的數(shù)值頻散進行了對比研究． 理論分析和數(shù)值計算結果表明： 有限元三角網(wǎng)格中波動的數(shù)值頻散受網(wǎng)格布局、 波傳播方向、 單元網(wǎng)格縱橫比以及質量矩陣的影響； 一致質量矩陣的數(shù)值頻散比集中質量矩陣更易受到波傳播方向的影響； 不合理的三角網(wǎng)格單元會對數(shù)值相速度（數(shù)值頻散）產生不良影響； 正三角網(wǎng)格中波動的數(shù)值頻散幾乎不受波傳播方向的影響； 一致質量矩陣與集中質量矩陣的線性組合能夠有效地壓制數(shù)值頻散．

有限元法 三角網(wǎng)格 波動方程 數(shù)值頻散 組合質量矩陣 數(shù)值仿真

引言

隨著高性能計算機的不斷發(fā)展， 有限元法在分析波的傳播問題上得到廣泛應用． 利用有限元方法模擬連續(xù)介質中的波動問題時， 需要將連續(xù)介質模型離散為網(wǎng)格模型， 這必然會引起精度上的誤差， 使得具有不同頻率的波動表現(xiàn)出不同的相速度（數(shù)值相速度）， 導致波場中出現(xiàn)偽波動， 即數(shù)值頻散（董良國， 李培明， 2004； 吳國忱， 王華忠， 2005）． 因此， 網(wǎng)格離散引起的數(shù)值頻散是影響波動方程數(shù)值計算精度的主要因素． 研究影響數(shù)值頻散的因素， 對提高波動方程數(shù)值模擬的精度具有指導意義．

波動方程有限元數(shù)值解的頻散分析要比有限差分法數(shù)值解的頻散分析難得多（Marfurt， 1984； 李勝軍等， 2008）， 這主要是因為有限元法具有多種網(wǎng)格剖分方式、 多種插值函數(shù)以及不同特性的質量矩陣（廖振鵬， 劉晶波， 1992； 薛東川等， 2007）， 這些因素都會對數(shù)值頻散產生影響． 波動方程有限元數(shù)值解的頻散問題一直是眾多研究者研究的熱點． 例如： Mullen和Belytschko（1982）最早對該問題進行研究， 考慮矩形網(wǎng)格以及不同形狀的三角網(wǎng)格， 得出了矩形網(wǎng)格比三角網(wǎng)格具有更高精度的結論； 宗福開（1984）分析波動方程有限元數(shù)值解的頻散特性， 給出了波動問題有限元離散化的準則， 并以幾個典型實例來說明有限元法所能達到的精確程度； Abboud和Pinsky（1992）在三維網(wǎng)格單元下， 對聲波方程有限元解的頻散進行研究， 討論了網(wǎng)格剖分及質量矩陣對數(shù)值頻散的影響； 房營光和莫海鴻（2000）采用含有頻率的高階位移函數(shù)對有限元網(wǎng)格中波動的頻散與穩(wěn)定性條件進行了改進； 薛東川和王尚旭（2008）從數(shù)值模擬角度研究了一致質量矩陣和集中質量矩陣對數(shù)值頻散的影響， 并采用二者的線性組合來壓制數(shù)值頻散； 此外， 波動方程譜元法數(shù)值解的頻散特性（De Basabe， Sen， 2007； Seriani， Oliveira， 2008； Liuetal， 2012）以及間斷有限元法數(shù)值解的頻散特性（Huetal， 1999； Ainsworthetal， 2006； De Basabeetal， 2008； 賀茜君等， 2014）也得到了廣泛研究．

本文針對有限元三角網(wǎng)格中波動的數(shù)值頻散問題， 推導了集中質量矩陣和一致質量矩陣的頻散函數(shù)， 通過頻散曲線對比分析了不同質量矩陣的數(shù)值頻散特性， 最后通過數(shù)值仿真對相關結論進行了驗證．

1 有限元三角網(wǎng)格單元的波動方程

使用伽勒金法求解二維標量波動方程得到有限元方程組（杜世通， 1982）為

（1）

式中，Ke和Me分別為單元剛度矩陣和單元質量矩陣，u為單元節(jié)點位移列向量． 計算出每個單元的剛度矩陣和質量矩陣后， 集成為結構剛度矩陣K和結構質量矩陣M， 得到結構的有限元方程組為

（2）

式中U為結構節(jié)點位移列向量．

圖1 三角網(wǎng)格單元示意圖

圖1給出了三角網(wǎng)格單元的示意圖. 三角網(wǎng)格單元的3個頂點按逆時針順序分別記為1， 2， 3， 其面積坐標可作如下定義（徐世浙， 1994）：

（3）

式中，a1=z3-z2，a2=z1-z3，a3=z2-z1，b1=x2-x3，b2=x3-x1，b3=x1-x2，c1=x3z2-x2z3，c2=x1z3-x3z1，c3=x2z1-x1z2，Δ=（a2b1-a1b2）/2為三角網(wǎng)格單元的面積． 可見，a1，a2，a3,b1,b2，b3,c1,c2，c3和Δ為只與單元節(jié)點坐標有關的常數(shù)， 因此L1，L2和L3是x和z的線性函數(shù)． 形函數(shù)矩陣N為

N=[L1L2L3]T,

（4）

則

（5）

三角網(wǎng)格單元的剛度矩陣Ke和一致質量矩陣Me分別為

（6）

（7）

（8）

對于式（2）， 時間因子采用二階中心差分格式進行求解得到

（9）

式中Δt為時間采樣間隔．

2 有限元三角網(wǎng)格中波動的數(shù)值頻散函數(shù)

考慮圖2a和圖2b所示的網(wǎng)格結構， 假定區(qū)域為均勻、 無邊界的， 且不考慮震源的加載． 圖2a是在矩形單元上添加對角線形成的三角單元； 圖2b則是在兩條平行線之間以等腰三角單元進行剖分， 底邊與x軸方向平行． 為了便于敘述， 將圖2a和圖2b所示的網(wǎng)格結構分別記為I類和Ⅱ類三角網(wǎng)格結構． 兩種網(wǎng)格結構在x軸和z軸方向的采樣間隔分別為Δx和Δz， 單元網(wǎng)格縱橫比為γ=Δz/Δx． 圖中帶下標的大寫字母為結構節(jié)點編號， 小寫字母為單元編號， 阿拉伯數(shù)字為單元內節(jié)點編號．

圖2 共點單元網(wǎng)格結構

2.1 Ⅰ類三角網(wǎng)格結構的頻散函數(shù)

（10）

根據(jù)剛度（質量）矩陣的組裝原理， 結構剛度矩陣K第J2行需要進行疊加計算的元素為

結構集中質量矩陣Ml第J2行的非零元素為

（12）

結構一致質量矩陣M第J2行的非零元素為

（13）

記節(jié)點J2的坐標為（mΔx，nΔz）， 根據(jù)平面波理論得到該點處的位移值um， n（孫成禹， 2007）為

um, n=Aexp[ik（mΔxcosθ+nΔzsinθ）-iωt],

（14）

式中，A為振幅，k為波數(shù)，ω為角頻率，θ為波傳播方向與x軸的夾角． 網(wǎng)格離散后的相速度cp為

（15）

對式（2）進行第J2行運算， 可以得到節(jié)點J2的位移與周圍相關節(jié)點位移的關系．

2.1.1 集中質量矩陣的頻散函數(shù)

在集中質量矩陣下， 對式（2）進行第J2行運算可以得到

um, nK（J2,J2）+um+1, nK（J2,K2）+um, n+1K（J2,J3）+um+1, n+1K（J2,K3）=0.

（16）

將式（11）、 （12）代入式（16）得到

（17）

再將式（14）代入式（17）， 整理得到

（18）

其中

（19）

在后續(xù)的推導中還會用到如下表達式：

（20）

最后將式（15）代入式（18）， 得到集中質量矩陣的頻散函數(shù)為

（21）

2.1.2 一致質量矩陣的頻散函數(shù)

同理對式（2）進行第J2行運算得到

2γ2（um-1, n+um+1, n）]=0.

（22）

仿照式（21）的推導過程， 可以得到一致質量矩陣的頻散函數(shù)為

（23）

2.2 Ⅱ類三角網(wǎng)格結構的頻散函數(shù)

（24）

根據(jù)剛度（質量）矩陣的組裝原理， 結構剛度矩陣K第J2行需要進行疊加計算的元素為

（25）

結構集中質量矩陣Ml第J2行的非零元素為

（26）

結構一致質量矩陣M第J2行的非零元素為

（27）

2.2.1 集中質量矩陣的頻散函數(shù)

對式（2）進行第J2行運算可以得到

（28）

仿照式（21）的推導過程， 得到集中質量矩陣的頻散函數(shù)為

（29）

2.2.2 一致質量矩陣的頻散函數(shù)

同理對式（2）進行第J2行運算可以得到

（30）

仿照式（21）的推導過程， 得到一致質量矩陣的頻散函數(shù)為

（31）

3 數(shù)值頻散特性分析

根據(jù)奈奎斯特（Nyquist）頻率可知， 當一個波長內空間采樣點少于2個時會出現(xiàn)假頻現(xiàn)象， 因此進行以下數(shù)值頻散分析時，kΔx的取值范圍為 0—π（孫成禹等， 2009； Liu， Sen， 2009）．

3.1 Ⅰ類三角網(wǎng)格結構中的頻散特性分析

當θ取不同值（γ=1）時，cp/c與kΔx的關系如圖3所示． 可以看出：

圖3 θ取不同值時cp/c隨kΔx的變化關系

1） 集中質量矩陣使相速度滯后于真實速度， 一致質量矩陣使相速度超前于真實速度， 二者對相速度起著相反的作用．

2） 隨著傳播角θ的增加， 一致質量矩陣的相速度比集中質量矩陣的相速度分散得更開， 即一致質量矩陣的數(shù)值頻散比集中質量矩陣更易受到波傳播方向的影響．

3） 無論是集中質量矩陣還是一致質量矩陣， 當kΔx的值接近于0時（空間采樣間隔足夠小），cp/c值幾乎為1， 此時數(shù)值頻散能夠得到較好地壓制， 且?guī)缀醪皇懿▊鞑シ较虻挠绊懀?/p>

4） 如果認為當所有傳播方向上的頻散誤差R=|1-cp/c|×100% ≤0.5%時， 可以忽略數(shù)值頻散的影響， 則集中質量矩陣滿足該不等式的條件是kΔx≤0.3481， 即單元網(wǎng)格的最大邊長Δxmin應該小于或等于λ/20（λ為波長）； 一致質量矩陣下滿足該不等式的條件是kΔx≤0.2291， 即Δxmin≤λ/27． 因此， 對空間采樣間隔而言， 集中質量矩陣壓制數(shù)值頻散的效率比一致質量矩陣的高．

下面研究單元網(wǎng)格縱橫比γ和波傳播方向θ對數(shù)值頻散的影響． 假定kΔx=π/2， 圖4給出了當γ取不同值時cp/c與θ的關系． 可以看出：

1） 在空間采樣間隔不是足夠小的情況下， 數(shù)值頻散因波傳播方向的不同而發(fā)生顯著的變化， 且這種變化是周期性的． 在集中質量矩陣下， 當γ=1時， 周期為π/2， 以π/4+kπ/2（k=0， 1， 2， 3, …）為對稱軸； 當γ≠1時， 相應的周期為π， 對稱軸為π/2+kπ（k=0， 1， 2， 3, …）． 在一致質量矩陣下， 周期均為π， 在周期內不具有對稱性．

圖4 γ取不同值時cp/c隨θ的變化關系

2） 當θ較小時， 集中質量矩陣的數(shù)值頻散幾乎不受單元網(wǎng)格縱橫比的影響， 而一致質量矩陣的數(shù)值頻散則不具有此特性．

3） 在集中質量矩陣下（γ=1）， 當波沿水平方向（θ=0°）和垂直方向（θ=90°）傳播時， 其數(shù)值頻散最嚴重； 當波沿θ=45°方向傳播時， 其數(shù)值頻散最小（圖4a中的γ=1.0曲線）， 且水平方向的數(shù)值頻散不受單元網(wǎng)格縱橫比的影響， 即通過減小單元網(wǎng)格縱橫比來壓制數(shù)值頻散不是一種行之有效的方法．

4） 在一致質量矩陣下， 水平方向的數(shù)值頻散也不受單元網(wǎng)格縱橫比的影響．

3.2 Ⅱ類三角網(wǎng)格結構中的頻散特性分析

仍然假定kΔx=π/2， 其數(shù)值頻散cp/c與θ及γ的關系分別如圖5和圖6所示． 可以看出：

1） 在集中質量矩陣下， 當單元網(wǎng)格縱橫比γ<0.5時， 得到了相速度有可能大于介質真實速度的結果， 該結果與集中質量矩陣使相速度滯后于真實速度相反． 進一步分析可以發(fā)現(xiàn)， 當γ=0.5時， 三角網(wǎng)格單元為等腰直角三角形； 當γ<0.5時， 三角網(wǎng)格單元為等腰鈍角三角形． 因此不難發(fā)現(xiàn)在集中質量矩陣下， 出現(xiàn)相速度大于介質真實速度的原因是單元剖分不合理， 出現(xiàn)了鈍角三角單元， 影響數(shù)值計算的精度． 同理， 當單元剖分不合理時， 一致質量矩陣的相速度突變得非常大， 對數(shù)值頻散會產生不良影響． 因此， 對這類網(wǎng)格結構而言， 選擇合適的網(wǎng)格縱橫比至關重要， 一定要避免出現(xiàn)過銳或過鈍的三角單元．

圖5 γ取不同值時cp/c隨θ的變化關系. （a） 集中質量矩陣； （b） 一致質量矩陣

圖6 θ取不同值時cp/c隨γ的變化關系. （a） 集中質量矩陣； （b） 一致質量矩陣

2） 在該類網(wǎng)格結構下， 兩種質量矩陣的數(shù)值頻散具有相同的周期性及對稱性．

3） 當單元網(wǎng)格縱橫比為0.866時， 不同傳播方向上的數(shù)值頻散幾乎相等（圖6中所有曲線的交點）， 此時三角網(wǎng)格單元為正三角形．

圖7中的γ=0.866曲線和圖8a分別展示了當單元為正三角形時數(shù)值頻散與波傳播方向及空間采樣間隔的關系． 可以看出， 圖7中的γ=0.866曲線周期為π/3， 即正三角網(wǎng)格下的周期由π變成了π/3， 因為正三角網(wǎng)格下的結構只需旋轉π/3就不會發(fā)生變化． 圖8a中的曲線基本重合， 也進一步證實了正三角單元下的數(shù)值頻散幾乎不受波傳播方向的影響. 圖8b則為單元網(wǎng)格縱橫比為0.4時數(shù)值頻散與空間采樣間隔的關系． 可以看出， 在集中質量矩陣下， 當波沿水平方向及其附近傳播時， 隨著空間采樣間隔的增加， 出現(xiàn)了相速度大于介質真實速度的現(xiàn)象， 且一致質量矩陣的相速度變得非常大．

圖7 γ取不同值時cp/c隨θ的變化關系. （a） 集中質量矩陣； （b） 一致質量矩陣

圖8 θ取不同值時cp/c隨kΔx的變化關系

3.3 一致質量矩陣與集中質量矩陣的線性組合對數(shù)值頻散的影響

（32）

式中α為組合系數(shù)， 則Ⅰ類網(wǎng)格結構的組合質量矩陣頻散函數(shù)可以表示為

（33）

Ⅱ類網(wǎng)格結構的組合質量矩陣頻散函數(shù)可以表示為

（34）

當組合系數(shù)α=1時， 組合質量矩陣的頻散函數(shù)變?yōu)橐恢沦|量矩陣的頻散函數(shù)； 當α=0時， 組合質量矩陣的頻散函數(shù)變?yōu)榧匈|量矩陣的頻散函數(shù)． 對于Ⅰ類網(wǎng)格結構， 當組合系數(shù)α=0.5、 單元網(wǎng)格縱橫比γ=1時， 組合質量矩陣的數(shù)值頻散與空間采樣間隔的關系如圖9a所示； Ⅱ類網(wǎng)格結構的組合系數(shù)α也取為0.5、 單元網(wǎng)格為正三角形時， 組合質量矩陣的數(shù)值頻散與空間采樣間隔的關系如圖9b所示．

圖9 θ取不同值時cp/c隨kΔx的變化關系

Ⅰ類網(wǎng)格結構采用組合質量矩陣后， 整體上數(shù)值頻散得到有效地壓制， 但某些傳播方向上的數(shù)值頻散也會增加． 例如圖9a中π/4方向上的數(shù)值頻散， 當kΔx=2時， 頻散誤差為10.1%， 而采用集中質量矩陣時， 相應的頻散誤差為8.1%． 此外， 還應注意到當采用組合質量矩陣后， 在某些傳播方向上相速度小于介質的真實速度， 而在另一些傳播方向上， 相速度則大于介質的真實速度， 這意味著相速度失去向真實速度收斂的單調性．

對于Ⅱ類網(wǎng)格結構（γ=0.866）， 采用組合質量矩陣后，cp/c接近于1的范圍得到顯著地延伸（圖9b）． 同樣認為當所有傳播方向上的頻散誤差R≤0.5%時， 可以忽略數(shù)值頻散的影響， 則在組合質量矩陣下kΔx≤1.421， 即一個波長內所包含的空間采樣點不應少于5個； 而單獨采用集中質量矩陣或一致質量矩陣時， 相應的采樣點不應少于15個． 最后取kΔx=π/2， 兩種網(wǎng)格結構的部分數(shù)值頻散如表1所示．

表1 kΔx=π/2時的頻散cp/cTable 1 Dispersion cp/c at kΔx=π/2

4 數(shù)值仿真

由式（9）可知， 用有限元法解波的傳播問題時， 需要在時間遞推過程中求解大型稀疏方程組． 結構一致質量矩陣M是大型稀疏矩陣， 對其求逆不僅需要占用大量機時， 而且引入數(shù)值計算誤差的可能性也增大； 而結構集中質量矩陣Ml是一對角矩陣， 在求解式（9）時， 也就避免了大型稀疏矩陣的求逆運算． 因此， 本文的數(shù)值仿真采用集中質量矩陣． 考慮縱波波速c=2.0 km/s的均勻介質模型， 計算區(qū)域為1.8 km×1.8 km， 震源采用主頻為33 Hz的雷克子波， 位于模型中央． 兩種網(wǎng)格結構的空間橫向采樣間隔Δx均為6 m， 時間采樣間隔Δt=0.5 ms， 改變單元網(wǎng)格縱橫比γ， 得到0.43 s時刻的波場快照如圖10所示． 圖10a （γ=1）和圖10b （γ=0.3）為在Ⅰ類網(wǎng)格結構下得到的波場快照． 可以看出: 當單元網(wǎng)格縱橫比γ=1時， 數(shù)值頻散在θ=0°時最大， 且隨著θ的增大而減小， 在θ=45°時最小， 這與圖4a中的γ=1.0曲線一致； 將單元網(wǎng)格縱橫比減小到0.3時， 只有當θ≥45°時， 數(shù)值頻散才得到有效地壓制． 在Ⅱ類網(wǎng)格結構下得到的波場快照分別如圖10c （γ=1）、 圖10d （γ=0.866）、 圖10e （γ=0.5）和圖10f （γ=0.3）所示． 可以看出： 當單元網(wǎng)格縱橫比由1減小到0.866時， 數(shù)值頻散在一定程度上得到減弱， 且此時各個傳播方向上的數(shù)值頻散基本一致； 當單元網(wǎng)格縱橫比減小到0.5時， 頻散壓制效果尤為明顯； 當單元網(wǎng)格縱橫比減小到0.3時， 在水平方向附近出現(xiàn)相速度大于介質真實速度的現(xiàn)象． 以上數(shù)值仿真結果均與理論分析一致．

圖10 不同條件下的數(shù)值頻散現(xiàn)象（0.43 s）

Fig.10 The numerical dispersion under different conditions （0.43 s）

（a） Ⅰ類網(wǎng)格結構（γ=1）; （b） Ⅰ類網(wǎng)格結構（γ=0.3）; （c） Ⅱ類網(wǎng)格結構（γ=1）;（d） Ⅱ類網(wǎng)格結構（γ=0.866）; （e） Ⅱ類網(wǎng)格結構（γ=0.5）; （f） Ⅱ類網(wǎng)格結構（γ=0.3）

（a） Mesh structure Ⅰ（γ=1）; （b） Mesh structure Ⅰ（γ=0.3）; （c） Mesh structure Ⅱ（γ=1）;（d） Mesh structure Ⅱ（γ=0.866）; （e） Mesh structure Ⅱ（γ=0.5）; （f） Mesh structure Ⅱ（γ=0.3）

5 結論

本文研究了有限元三角網(wǎng)格中不同質量矩陣下波動的數(shù)值頻散特性， 得到如下結論：

1） 在非鈍角三角網(wǎng)格中， 集中質量矩陣使數(shù)值相速度滯后于真實速度， 一致質量矩陣使數(shù)值相速度超前于真實速度．

2） 一致質量矩陣的數(shù)值頻散比集中質量矩陣更易受到波傳播方向的影響．

3） 在Ⅰ類三角網(wǎng)格結構中， 水平方向（θ=0°）的數(shù)值頻散不受單元網(wǎng)格縱橫比的影響， 集中質量矩陣與一致質量矩陣的線性組合可以減弱大部分傳播方向上的數(shù)值頻散， 但在某些傳播方向上， 數(shù)值頻散反而增加； 此外， 采用組合質量矩陣后， 因傳播方向的不同， 相速度有可能小于介質的真實速度， 也有可能大于介質的真實速度．

4） 在Ⅱ類三角網(wǎng)格結構中， 過銳或過鈍的三角單元不僅會導致集中質量矩陣的相速度大于介質的真實速度， 而且還會使一致質量矩陣的相速度變得異常的大， 對數(shù)值解的精度產生不良影響， 甚至會導致錯誤的結果， 而正三角單元中的數(shù)值頻散受波傳播方向的影響較?。?因此在網(wǎng)格剖分中， 單元的形狀接近等邊三角形最好， 不宜過銳或過鈍． 在保證三角單元非過銳或過鈍的前提下， 減小單元網(wǎng)格縱橫比可以顯著地壓制各個傳播方向上的數(shù)值頻散．

5） 在正三角網(wǎng)格結構中采用組合質量矩陣后， 有效壓制數(shù)值頻散的最大空間采樣間隔可以得到顯著增加， 且不會出現(xiàn)相速度在某些傳播方向上小于介質的真實速度， 而在另一些傳播方向上大于介質真實速度的現(xiàn)象．

董良國， 李培明. 2004. 地震波傳播數(shù)值模擬中的頻散問題[J]. 天然氣工業(yè)， 24（6）: 53--56.

Dong L G， Li P M. 2004. Dispersive problem in seismic wave propagation numerical modeling[J].NaturalGasIndustry， 24（6）: 53--56 （in Chinese）.

杜世通. 1982. 變速不均勻介質中波動方程的有限元法數(shù)值解[J]. 華東石油學院學報， 6（2）: 1--20.

Du S T. 1982. Finite element numerical solution of wave propagation in non-homogeneous medium with variable velocities[J].JournalofEastChinaPetroleumInstitute， 6（2）: 1--20 （in Chinese）.

房營光， 莫海鴻. 2000. 有限元網(wǎng)格中波動的頻散與穩(wěn)定性的一種改進方法[J]. 地震工程與工程振動， 20（1）: 21--26.

Fang Y G， Mo H H. 2000. An improved method for dispersion and stability of wave motion in finite element meshes[J].EarthquakeEngineeringandEngineeringVibration， 20（1）: 21--26 （in Chinese）.

賀茜君， 楊頂輝， 吳昊. 2014. 間斷有限元方法的數(shù)值頻散分析及其波場模擬[J]. 地球物理學報， 57（3）: 906--917.

He X J， Yang D H， Wu H. 2014. Numerical dispersion and wave-field simulation of the Runge-Kutta discontinuous Galerkin method[J].ChineseJournalofGeophysics， 57（3）: 906--917 （in Chinese）.

李勝軍， 孫成禹， 高建虎， 劉軍迎. 2008. 地震波數(shù)值模擬中的頻散壓制方法分析[J]. 石油物探， 47（5）: 444--448.

Li S J， Sun C Y， Gao J H， Liu J Y. 2008. Analysis of dispersion suppression in wave equation numerical simulation[J].GeophysicalProspectingforPetroleum， 47（5）: 444--448 （in Chinese）.

廖振鵬， 劉晶波. 1992. 波動有限元模擬的基本問題[J]. 中國科學： B輯， 22（8）: 874--882.

Liao Z P， Liu J B. 1992. The basic problems of wave motion with finite element simulation[J].ScienceinChina:SeriesB， 22（8）: 874--882 （in Chinese）.

孫成禹. 2007. 地震波理論與方法[M]. 東營: 中國石油大學出版社: 31--37.

Sun C Y. 2007.TheoryandMethodsofSeismicWaves[M]. Dongying: China University of Petroleum Press: 31--37 （in Chinese）.

孫成禹， 宮同舉， 張玉亮， 張文穎. 2009. 波動方程有限差分法中的頻散與假頻分析[J]. 石油地球物理勘探， 44（1）: 43--48.

Sun C Y， Gong T J， Zhang Y L， Zhang W Y. 2009. Analysis on dispersion and alias in finite-difference solution of wave equation[J].OilGeophysicalProspecting， 44（1）: 43--48 （in Chinese）.

吳國忱， 王華忠. 2005. 波場模擬中的數(shù)值頻散分析與校正策略[J]. 地球物理學進展， 20（1）: 58--65.

Wu G C， Wang H Z. 2005. Analysis of numerical dispersion in wave-field simulation[J].ProgressinGeophysics， 20（1）: 58--65 （in Chinese）.

徐世浙. 1994. 地球物理中的有限單元法[M]. 北京: 科學出版社: 39--42.

Xu S Z. 1994.FiniteElementMethodforGeophysics[M]. Beijing: Science Press: 39--42 （in Chinese）.

薛東川， 王尚旭， 焦淑靜. 2007. 起伏地表復雜介質波動方程有限元數(shù)值模擬方法[J]. 地球物理學進展， 22（2）: 522--529.

Xue D C， Wang S X， Jiao S J. 2007. Wave equation finite-element modeling including rugged topography and complicated medium[J].ProgressinGeophysics， 22（2）: 522--529 （in Chinese）.

薛東川， 王尚旭. 2008. 利用組合質量矩陣壓制數(shù)值頻散[J]. 石油地球物理勘探， 43（3）: 318--320.

Xue D C， Wang S X. 2008. Using combined mass matrix to suppress numeric dispersion[J].OilGeophysicalProspecting， 43（3）: 318--320 （in Chinese）.

宗福開. 1984. 波傳播問題中有限元分析的頻散特性及離散化準則[J]. 爆炸與沖擊， 4（4）: 16--23.

Zong F K. 1984. Frequency-dispersion characteristics and discretization of the finite element analysis in wave propagation problems[J].ExplosionandShockWaves， 4（4）: 16--23 （in Chinese）.

Abboud N N， Pinsky P M. 1992. Finite element dispersion analysis for the three-dimensional second-order scalar wave equation[J].IntJNumerMethodsEng， 35（6）: 1183--1218.

Ainsworth M， Monk P， Muniz W. 2006. Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation[J].JSciComput， 27（1/2/3）: 5--40.

De Basabe J D， Sen M K. 2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations[J].Geophysics， 72（6）: T81--T95.

De Basabe J D， Sen M K， Wheeler M F. 2008. The interior penalty discontinuous Galerkin method for elastic wave pro-pagation: Grid dispersion[J].GeophysJInt， 175（1）: 83--93.

Hu F Q， Hussaini M Y， Rasetarinera P. 1999. An analysis of the discontinuous Galerkin method for wave propagation problems[J].JComputPhys， 151（2）: 921--946.

Liu T， De Basabe J D， Li L， Hu T Y， Sen M K， Wei X C. 2012. Grid dispersion and stability of the spectral element method with triangular elements[C]∥82thAnnualInternationalMeeting，SEG，TechnicalProgramExpandedAbstracts. Las Vegas： Society of Exploration Geophysicists: 1--5.

Liu Y， Sen M K. 2009. A new time-space domain high-order finite-difference method for the acoustic wave equation[J].JComputPhys， 228（23）: 8779--8806.

Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J].Geophysics， 49（5）: 533--549.

Mullen R， Belytschko T. 1982. Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation[J].IntJNumerMethodsEng， 18（1）: 11--29.

Seriani G， Oliveira S P. 2008. Dispersion analysis of spectral element methods for elastic wave propagation[J].WaveMotion， 45（6）: 729--744.

Dispersion analysis and numerical simulation of wave motion in finite element algorithm with triangular meshes

1）SchoolofGeosciences，ChinaUniversityofPetroleum，ShandongQingdao266580，China2）WesternNewProspectResearchInstitute，ShengliOilfieldCompany，SINOPEC，ShandongDongying257061，China3）SchoolofGeodesyandGeomatics，WuhanUniversity，Wuhan430079，China

The finite element discretization of wave motion can cause numerical error. The numerical dispersion is in essence the non-physical solution caused by the propagation of numerical error. Not only has the numerical dispersion no practical meaning， but also it can affect the understanding of real fluctuations. In order to clarify the influence factors of numerical dispersion in finite element algorithm with triangular meshes， the dispersion functions of lumped mass matrix and consistent mass matrix are derived respectively， and the dispersion function of combined mass matrix is also given. And then the numerical dispersion of different mass matrices are compared. The results of theoretical analysis and numerical simulation indicate: ① The numerical dispersion in finite element algorithm with triangular meshes depend on the mesh layout， direction of wave propagation， the ratio of vertical length to horizontal length and the mass matrix; ② The numerical dispersion of consistent mass marix is more easily affected by wave propagation direction than that of lumped mass matrix; ③ The irrational triangular meshes have a bad effect on numerical phase velocity （numerical dispersion）; ④ Equilateral triangular meshes can provide superior results regardless of the propagation direction; ⑤ The linear combination of lumped mass matrix and consistent mass matrix can effectively suppress numerical dispersion.

finite element method; triangular mesh; wave equation; numerical dispersion; combined mass matrix; numerical simulation

10.11939/jass.2015.03.012.

國家“十二五”重點科技攻關項目（2011ZX05002-002）與中石化科技重大攻關項目（P13020）聯(lián)合資助.

2014-05-25收到初稿， 2014-08-08決定采用修改稿.

e-mail: pillar1979@163.com

10.11939/jass.2015.03.012

P315.3+1

A

王圣柱, 周建科, 張奎華, 王藝豪. 2015. 有限元三角網(wǎng)格中波動的頻散分析及數(shù)值仿真. 地震學報, 37（3）: 493--507.

Wang S Z, Zhou J K, Zhang K H, Wang Y H. 2015. Dispersion analysis and numerical simulation of wave motion in finite element algorithm with triangular meshes.ActaSeismologicaSinica, 37（3）: 493--507. doi:10.11939/jass.2015.03.012.