楊在林， 孫 鋮， 王 耀， 李志東

(哈爾濱工程大學航天與建筑工程學院，黑龍江 哈爾濱 150001)

含方形凹陷半無限非均勻介質波動問題FDM模擬①

楊在林， 孫 鋮， 王 耀， 李志東

(哈爾濱工程大學航天與建筑工程學院，黑龍江 哈爾濱 150001)

采用規(guī)則網(wǎng)格有限差分方法對二維平面彈性波動方程進行差分離散，得到相應的彈性波動方程的有限差分方程，再將彈性波動方程的差分格式與吸收邊界、自由邊界的離散形式結合形成彈性波動方程有限差分方程解決問題的主體，將其應用于含方形凹陷半無限非均勻介質的模型中進行數(shù)值模擬，得到此離散化模型中不同時刻不同節(jié)點的位移值。針對具體算例，運用上述方法結合科學計算軟件MATLAB和結果后處理軟件DIFEM ISOLINE PLOTER得到不同時刻的水平方向位移等值線圖與接收器測量點處的合成位移記錄，討論非均勻介質、吸收邊界、方形凹陷等對波動特性的影響。

非均勻介質； 波動； 有限差分方法； 凹陷； Higdon吸收邊界

0 引言

波動在聲學、電磁學、光學、地球物理學等領域有著重要的理論意義與應用價值。例如在地球物理學領域中，地震勘探(地震波勘探)是一種在石油工業(yè)領域中應用極其廣泛的勘探手段，其原理是通過觀測和分析地下介質密度與彈性模量的不同形成對人工激勵的不同響應，然后通過數(shù)值模擬得到的數(shù)據(jù)合理地推斷出巖層的巖性、構造等相關信息[1]。隨著科學技術的進步，現(xiàn)階段地震波數(shù)值模擬主要有三類方法：積分方程法、幾何射線法、波動方程法。其中波動方程法可以通過波動方程模擬出豐富的波動信息，能更好地解釋復雜地層彈性波的傳播機理。因此波動方程數(shù)值模擬在地震波模擬領域具有重要價值[2]。

彈性波動問題的研究在數(shù)值方法方面已經發(fā)展出很多的解法[3-9]，其中有限差分方法是求解波動方程等雙曲型偏微分方程的最常用的數(shù)值方法之一，具有運算速度快、編程簡單等優(yōu)點。通過推導得到彈性波動方程的有限差分格式，以及初始條件及自由邊界、吸收邊界的離散形式[10-12]，然后通過有限差分方程組的求解得到彈性波動問題的數(shù)值解[13-14]。

1 理論推導

1.1 波動方程的有限差分離散

假設彈性波在x-z平面內傳播，則其參數(shù)與y坐標無關，取正應力分量σx=τxx，σz=τzz，剪應力分量為εx=εxx，εz=εzz，剪應變分量為εxz。

平面問題的彈性動力方程為

(1)

幾何方程為

(2)

本構方程為

(3)

上式中

(4)

彈性系數(shù)矩陣為

(5)

則

(6)

將式(6)代入式(1)并展開得

(7)

基于規(guī)則網(wǎng)格的八節(jié)點顯式差分格式(圖1)。

圖1 八節(jié)點規(guī)則網(wǎng)格示意圖Fig.1 Diagram of regular grid with eight-node

令x方向上位移為μ，z方向上位移為ω，可以得到彈性波動方程的有限差分格式：

式(8)中

(9)

(10)

(11)

(12)

將上面推導出的差分公式和初始-邊值離散條件結合，可逐層推算出整個求解域的位移波場值，進而能夠進一步得到擁有復雜界面的非均勻各向異性介質彈性波場。

出平面波動方程形式為

(13)

可得到出平面波的格式：

(14)

式中uy(i，j，n)表示t=n時刻的坐標x=i，z=j處的節(jié)點位移。通過與初始-邊值條件聯(lián)立，可以求得各個時刻各個節(jié)點的出平面波的波幅。

1.2 自由邊界和吸收邊界

以水平方向上的自由邊界為例，為了求解自由邊界上的相關位移，可以在自由邊界上面再附加一層虛擬的邊界點，通過相關式子用已知參數(shù)代替虛擬點的參數(shù)，再帶入彈性波動方程有限差分方程中，從而求得自由邊界上的相關參數(shù)。由自由邊界上具有垂直于自由邊界的應力為零這一特性，可知：

(15)

通過中心差分方法，將應力方程離散化，得到離散方程：

(16)

對于水平的自由邊界，7點的相關位移可以用其他點代替。同理，8、6點的位移也可以表示出來，三點的差分格式為：

(17)

吸收邊界采用Higdon人工吸收邊界條件，其二階吸收邊界差分算子為：

(18)

對其進行差分離散、歸納，可以得到吸收邊界上的波場值表達式。彈性波動方程有限差分方程、自由邊界、吸收邊界與初始-邊值條件結合，從而能夠對問題進行求解。

2 穩(wěn)定性分析

基于二維均勻介質中的顯式格式的譜分析方法得到穩(wěn)定性條件，可以表示為：

(19)

式中υP表示縱波波速；Δx、Δz分別表示x、z方向上的空間步長。由式(19)可知，穩(wěn)定性條件與橫波波速無關。通過數(shù)學方法可將此關系推廣到非均勻質的規(guī)則網(wǎng)格差分格式中，其更一般形式為：

(20)

3 算例分析與討論

3.1 模型的建立

設模型尺寸1 500 m×500 m，震源位于(490 m，190 m)處，震源函數(shù)為

(21)

圖2 模型網(wǎng)格劃分示意圖Fig.2 Mesh of the model

表1 模型規(guī)則網(wǎng)格差分與非均勻介質相關系數(shù)

3.2 波場位移的等值線圖

由不同時刻的水平方向位移等值線圖(圖3)可以發(fā)現(xiàn)，波傳播到方形凹陷時將發(fā)生反射與散射。在非均勻介質中波在不同位置的傳播速度不同，因此在震源上下方向呈現(xiàn)出不對稱性；同時介質的彈性系數(shù)矩陣和密度不同，在距震源相同距離的上下位置處位移響應值也不同；波在吸收邊界處的吸收效果較理想，說明二階Higdon吸收邊界適用于非均勻介質波動問題的求解。

圖3 水平方向位移等值線圖Fig.3 Displacement contour map in horizontal direction

3.3 接收器測量點處的合成位移記錄

通過接收器記錄的數(shù)據(jù)(圖4)可以發(fā)現(xiàn)，以x方向位移為例，波在具有方形凹陷的介質里傳播時，方形凹陷能對凹陷后方波的幅值產生一定的削弱作用。這說明在建筑物的選址過程中，需要考慮周邊凹陷地形(例如溝壑、湖泊等)對建筑物抗震性能的影響。

4 結論

本文是基于彈性波動理論，通過有限差分數(shù)值方法將彈性波動方程離散化得到相應的有限差分計算格式，結合初始-邊值條件，求得各點在各個時刻的波場位移值，通過對計算結果進行分析，可以得到以下結論：

圖4 不同觀測點處接收器x、z方向位移記錄Fig.4 Displacement records in x and z direction at different points

(1) 通過觀察波在吸收邊界處的吸收效果，發(fā)現(xiàn)在適當?shù)臈l件參數(shù)下，二階Higdon吸收邊界在非均勻介質波動問題的求解方面具有較好的效果，說明該人工邊界能夠適用于非均勻介質波動問題。

(2) 波在具有方形凹陷的介質里傳播時，方形凹陷能對凹陷后方波的強度產生一定的削弱作用。這對優(yōu)化建筑選址、提高建筑物抗震能力具有一定的參考價值。

References)

[1] 孫衛(wèi)濤.彈性波動方程的有限差分數(shù)值方法[M].北京：清華大學出版社，2009.SUN Wei-tao.Finite Element Method of Equations for Elastic Waves[M].Beijing：Tsinghua University Press，2009.(in Chinese)

[2] 裴正林，牟永光.地震波傳播數(shù)值模擬[J].地球物理學進展，2004，19(4)：933-941.PEI Zheng-lin，MU Yong-guang.Numerical Simulation of Seismic Wave Propagation[J].Progress in Geophysics，2004，19(4)：933-941.(in Chinese)

[3] Alterman Z，Karal F C JR.Propagation of Elastic Waves in Layered Media by Finite Difference Methods[J].Bulletin of the Seismological Society of America，1968，58(1)：367-398.

[4] Akik K，Richards P G.Quantitative Seismology[M].New York：W H Freeman and Company，1980.

[5] Smith W D.The Application of Finite-element Analysis to Body Wave Propagation Problem[J].Geophysics J Roy Astr Soc，1975，42：747-768.

[6] Cerveny V，Molotkov I A，Psencik I.Ray Method in Seismology[M].Praha：Univ Karlova1977.

[7] 李信富，李小凡，張美根.地震波數(shù)值模擬方法研究綜述[J].防災減災工程學報，2007，27(2)：241-248.LI Xin-fu，LI Xiao-fan，ZHANG Mei-gen.Review of Seismic Wave Numerical Modeling Methods[J].Journal of Disaster Prevention and Mitigation Engineering，2007，27(2)：241-248.(in Chinese)

[8] DUAN Yu-ting，HU Tian-yue，YAO Feng-chang，et al.3D Elastic Wave Equation Forward Modeling Based on the Precise Integration Method[J].Applied Geophysics，2013：71-78，118-119.

[9] 廖振鵬.工程波動理論導論(第二版)[M].北京：科學出版社，2002.LIAO Zhen-peng.Introduction to Engineering Wave Theory(the Second Edition)[M].Beijing：Science Press，2002.(in Chinese)

[10] Higdon R L.Absorbing Boundary Conditions for Difference Approximations to The Multidimensional Wave Equation[J].Math Comp，1986，47(176)：437-459.

[11] Higdon R L.Numerical Absorbing Boundary Conditions for the Wave Equation[J].Math Comp，1987，49(179)：65-90.

[12] 任濟時.高階Higdon吸收邊界條件的直接算法及其評估[J].電子學報，1997，25(5)：110-112.REN Ji-shi.A Directly Implemented Method and Evaluation for the Higher Order Higdon Absorbing Boundary Condition[J].Acta Electronica Sinica，1997，25(5)：110-112.(in Chinese)

[13] 孫衛(wèi)濤，楊慧珠，舒繼武.非均勻介質彈性波動方程的不規(guī)則網(wǎng)格有限差分方法[J].計算力學學報，2004，21(2)：135-142.SUN Wei-tao，YANG Hui-zhu，SHU Ji-wu.Finite Difference Method of Lrregular Grid for Elastic Wave Equation in Heterogeneous Media[J].Chinese Journal of Computational Mechanics，2004，21(2)：135-142.(in Chinese)

[14] 唐文,王尚旭,袁三一.起伏地表二階彈性波方程差分策略穩(wěn)定性分析[J].石油物探，2013，52(5)：457-463，441.TANG Wen，WANG Shang-xu，YUAN San-yi.Stability Analysis of Differential Strategy for Rugged Topography by Second-order Elastic Wave Equation Based on Coordinate Transformation[J].Geophysical Prospecting for Petroleum，2013，52(5)：457-463，441.(in Chinese)

FDM Simulation of Wave Motion in a Semi-infinite Inhomogeneous Medium with a Rectangular Depression

YANG Zai-lin， SUN Cheng， WANG Yao， LI Zhi-dong

(CollegeofAerospaceandCivilEngineering，HarbinEngineeringUniversity，Harbin，Heilongjiang150001，China)

Using the finite difference method，the two-dimensional elastic wave equation is processed to obtain the finite difference equation.Then，the differential elastic wave equation format is combined with the absorbing boundary and the discrete form of the free boundary.The elastic wave equation finite difference equation is formed to solve the problem.Using the finite difference equation，the displacement on different nodes can be obtained at different times for the rectangle-shaped depression of a semi-infinite inhomogeneous medium.Using the finite difference equation of the two-dimensional elastic wave equation combined with the scientific computation software，MATLAB，and the post-processor software，DIFEM ISOLINE PLOTER，the displacement isoline of the rectangle-shaped depression of the semi-infinite inhomogeneous medium in thexdirection can be obtained for different times.The effect of the inhomogeneous media，absorbing boundary，and the rectangle-shaped depression on the wave characteristics is analyzed numerically.

inhomogeneous media； wave motion； finite difference method； depression； Higdon absorbing boundary

2014-08-20

2015年地震行業(yè)科研專項經費項目(201508026-02)；黑龍江省自然科學基金(A201310)；黑龍江省博士后科研啟動金(LBH-Q13040)

楊在林(1971-)，男，教授，博士生導師，主要從事彈性波動理論的相關研究.E-mail：yangzailin00@163.com

O343.7

A

1000-0844(2015)02-0553-06

10.3969/j.issn.1000-0844.2015.02.0553