李 志， 薛霆虓,2,張海亮

(1.桂林理工大學(xué) 地球科學(xué)學(xué)院，廣西 桂林 541004；2.桂林理工大學(xué) 廣西地質(zhì)工程中心區(qū)重點實驗室，廣西 桂林 541004；3.中化地質(zhì)礦山總局 山東地質(zhì)勘查院，山東 濟南 250013)

?

有限元地震數(shù)值模擬中無限元邊界的應(yīng)用

李 志1， 薛霆虓1,2,張海亮3

(1.桂林理工大學(xué) 地球科學(xué)學(xué)院，廣西 桂林 541004；2.桂林理工大學(xué) 廣西地質(zhì)工程中心區(qū)重點實驗室，廣西 桂林 541004；3.中化地質(zhì)礦山總局 山東地質(zhì)勘查院，山東 濟南 250013)

在用有限元做地震數(shù)值模擬中，對于研究無限區(qū)域的問題，需要在模型的邊界上采用一定的邊界條件，其他方法設(shè)置邊界條件較為困難，而利用無限元方法，可以比較有效地設(shè)計吸收邊界，編程也較為方便。本文設(shè)計了一些簡單的模型，利用有限元方法進行模擬，在邊界處設(shè)計了無限元，并得到了地震快照和人工地震記錄。分析對比結(jié)果可以明顯地看出，在使用無限元的方法后模型邊界的影響幾乎被消除，且適用于復(fù)雜地形。在解決無限域的問題上使用無限元能更真實地模擬地震波傳播過程。

地震數(shù)值模擬；有限元；無限元

1 引 言

地震數(shù)值模擬是地震勘探和地震學(xué)的重要基礎(chǔ)。所謂地震數(shù)值模擬就是對實際的地質(zhì)、 地球物理問題作適當?shù)慕疲纬珊喕臄?shù)學(xué)模型，采用數(shù)值計算方法獲取地震響應(yīng)的過程。地震數(shù)值模擬是了解地震波在地下介質(zhì)中傳播規(guī)律、幫助解釋觀測數(shù)據(jù)的有效手段[1]。這種地震數(shù)值模擬方法已在地震勘探和天然地震領(lǐng)域中得到廣泛的應(yīng)用。它不但在石油、天然氣、煤、金屬和非金屬等礦產(chǎn)資源及工程和環(huán)境地球物理中得到普遍的應(yīng)用，而且在地震災(zāi)害預(yù)測、地震區(qū)帶劃分以及地殼構(gòu)造和地球內(nèi)部結(jié)構(gòu)研究中，也得到相當廣泛的應(yīng)用[2]。

地震數(shù)值模擬有多種方法，比如有限差分法、有限元法等等。其中，有限差分法編程相對簡單，計算成本少，計算速度快，因此在地震勘探等領(lǐng)域得到了廣泛的應(yīng)用。但是有限差分法對復(fù)雜模型，比如地形起伏等情況，編程較為困難。有限元法雖然編程較為繁瑣，計算成本較大，但是可以模擬復(fù)雜模型，比如有起伏地形的模型。有限元法也是地震數(shù)值模擬中的一種較好的方法。然而，在使用有限元方法過程中，只能通過建立有限域的模型來模擬無限域的地下介質(zhì)，人為地引入介質(zhì)邊界，從而使得計算結(jié)果中出現(xiàn)無實際意義的邊界強反射，嚴重地干擾了有效波信息，并影響了對地震波傳播規(guī)律的正確認識[3]。為此，需要解決如何定義無限的區(qū)域，通常人為地截取足夠大的區(qū)域進行有限元網(wǎng)格劃分[4]，或者在模型的邊界上采用一定的邊界條件從而達到消減人為邊界反射的目的[5]。

在有限差分方法中，目前使用廣泛的兩類吸收邊界條件：海綿吸收邊界條件和傍軸近似吸收邊界條件。海綿吸收邊界條件利用黏滯邊界或靠近邊界的條帶范圍內(nèi)對入射波進行衰減[6,7]。Berenger[8]提出的理想匹配層法( PML method)是對海綿吸收邊界條件的進一步改善。傍軸近似吸收邊界條件基于單程波傳播的方法原理，將不同近似條件下的單程波方程作為不同邊界處的吸收邊界條件[9-12]。在有限元中使用PML等方法相對較為復(fù)雜，編程也困難。

在有限元的模型邊界中加入無限元，可以完美實現(xiàn)無限域的模擬且編程較為容易。1973年R Ungless[13]首先提出無限元的思想，此后Bettess[14]、Beer和Meek[15]、Zienkiewicz[16]等人發(fā)展了這種方法[17]。有限元和無限元耦合模型在求解實際問題方面有著廣泛的實用性，在模擬無限域問題上有很大的優(yōu)勢[19-23]。

2 無限元的原理

用有限元和無限元相結(jié)合的方法求解問題時，除了要用一般的有限元構(gòu)造近場的網(wǎng)格外，還需要用無限元構(gòu)造遠場網(wǎng)格。有限元與無限元之間存在分界點，通?？杉僭O(shè)該分界點到無限處的位移是線性分布的，這個線性假設(shè)相當合理，該無限元網(wǎng)格對近場網(wǎng)格的影響的精度是足夠的，為此我們需要構(gòu)造無限元的統(tǒng)一形式的基本解u。

以一維問題為例，如圖1中所示，設(shè)有限元和無限元的分界點為r1，奇異點位置為r0，r0到r1的距離為a，單位為m，r1作為分界點，位移為u1；r2為延伸到無限元中的一點，到分界點r1的距離為a，其位移為u2;r3作為無限遠點，位移自然為0；將坐標進行如下轉(zhuǎn)換。

圖1 無限元中的極點位置Fig.1 The pole position in infinite element

(1)

式中：s是映射坐標系的坐標；r是距離極點的距離。

在無限元中采用標準的二次位移插值函數(shù)：

(2)

或者寫成

(3)

這樣在1點，r=a,s=-1,u=u1；在2點，r=2a,s=0,u=u2；在3點，r=∞,s=1,u=0，說明構(gòu)造的位移插值函數(shù)滿足要求[19]。

3 算例分析

3.1 模型設(shè)計

本文選取的是二維彈性模型進行分析，有限元計算的區(qū)域范圍是40 m×20 m，材料的楊氏模量為3×1010Pa，泊松比為0.2，密度2 500 kg/m3。本文中的3個算例都是在模型的坐標原點處加載一個震源，震源模型為雷克子波，波形圖如圖2所示。震源作用時長為2 ms，作用在原點左右1 m的范圍上，振幅為500 000 Pa，如圖3所示。本文分別就有限元模型和有限元與無限元耦合模型討論在30 ms內(nèi)各節(jié)點的位移和速度情況。

本文的模型中使用了0.5 m的網(wǎng)格和0.1 m的加密網(wǎng)格(在震源作用半徑1 m的半圓中加密網(wǎng)格)。在有限元模型中共使用了4 511個四邊形網(wǎng)格和142個三角形網(wǎng)格，如圖4所示。在有限元和無限元耦合模型中共使用了4 581個四邊形網(wǎng)格、128個三角形網(wǎng)格和160個無限單元網(wǎng)格，如圖4、圖5所示。為了更好地說明無限元的實用性，筆者還建立了一個有地形起伏的模型，使用了4 681個四邊形網(wǎng)格、147個三角形網(wǎng)格和164個無限單元網(wǎng)格，如圖6所示。

圖2 雷克子波波形Fig.2 The chart of Ricker wavelet

圖3 加載震源Fig.3 The figure of loading source

圖4 模型的有限元網(wǎng)格Fig.4 The finite element mesh of the model

圖5 模型的有限元和無限元耦合網(wǎng)格Fig.5 The model mesh with finite element and infinite element

圖6 含地形的有限元和無限元耦合網(wǎng)格Fig.6 The terrain model mesh with finite element and infinite element

3.2 結(jié)果分析

對上述三種模型應(yīng)用有限元軟件做數(shù)值模擬計算，可以得到下面一系列圖形，圖7(a)～圖7

(e)分別表示有限元網(wǎng)格的模型在時間為2 ms、4 ms、6 ms、8 ms和10 ms時模型位移場的場值。圖8(a)～圖8(e)分別表示有限元和無限元耦合網(wǎng)格的模型在時間為2 ms、4 ms、6 ms、8 ms和10 ms時模型位移場的場值。圖9(a)～圖9(e)分別表示含地形的有限元和無限元耦合網(wǎng)格的模型在時間為2 ms、4 ms、6 ms、8 ms和10 ms時模型位移場的場值。

圖7(a)～圖7(e)、圖8(a)～圖8(e)和圖9(a)～圖9(e)分別表示地震波向外傳播的一個過程，但是由上述多幅圖可以看出，在地震波還未傳達到模型邊界的時候，即在6 ms以前，三種模型的波形圖完全一致。但在6 ms的時刻未使用無限元的模型，在地面明顯出現(xiàn)了反射波，8 ms和10 ms的波形圖更加明顯，而使用了無限元的模型完全看不到模型邊界所產(chǎn)生的反射波。對比圖8(a)～圖8(e)和圖9(a)～圖9(e)可以看出使用無限元對于含地形起伏的模型一樣地適用。

圖7 2 ms、4 ms、6 ms、8 ms和10 ms時位移場Fig.7 Displacement field in 2 ms, 4 ms, 6 ms, 8 ms and 10 ms

圖8 2 ms、4 ms、6 ms、8 ms和10 ms時位移場Fig.8 Displacement field in 2 ms, 4 ms, 6 ms, 8 ms and 10 ms

圖9 2 ms、4 ms、6 ms、8 ms和10 ms時位移場Fig.9 Displacement field in 2 ms, 4 ms, 6 ms, 8 ms and 10 ms

為了更加清楚地對比這三種模型的差異，在地表每隔2m選取一個節(jié)點垂直方向上的位移隨時間變化的數(shù)據(jù)，并運用GMT畫圖軟件，畫出20個節(jié)點的地震時間剖面圖，縱坐標表示時間長度共0.03 s，如圖10、圖11和圖12所示。

圖10 有限元網(wǎng)格的地震時間剖面Fig.10 Seismic time profile of finite element mesh

圖11 有限元與無限元耦合網(wǎng)格的地震時間剖面Fig.11 Seismic time profile of finite element and infinite element mesh

圖12 含地形的無限元與有限元耦合網(wǎng)格的地震時間剖面Fig.12 Seismic time profile of the terrain model of finite element and infinite element mesh

圖10與圖11和圖12對比而言可以清晰地看見底面的反射波，而圖11和圖12由于模型使用了無限元幾乎看不到反射波，而且圖12能看到地表地形起伏對初始波到來的時間的影響。

4 結(jié) 論

本文主要運用了有限元軟件做了一個地震數(shù)值模擬，對有限元模型和有限元和無限元耦合模型做了詳細的分析和討論，對比了這三種模型的結(jié)果，可以得到一個明顯的結(jié)論，在模擬無限域的問題上，使用無限元和有限元耦合模型可以很好地消除模型邊界的反射，更加接近于真實情況。有限元和無限元的結(jié)合，既能模擬復(fù)雜介質(zhì)，又能簡便地設(shè)置吸收邊界，達到模擬無限區(qū)域的目的。

[1]任志明,劉洋.一階彈性波方程數(shù)值模擬中的混合吸收邊界條件[J].地球物理學(xué)報,2014,57(2):595-606.

[2]牟永光,裴正林.三維復(fù)雜介質(zhì)地震數(shù)值模擬[M].北京：石油工業(yè)出版社，2005:1-2.

[3]馬在田.計算地球物理學(xué)概論[M].上海：同濟大學(xué)出版社,1997:10-71.

[4]戚玉亮,大塚久哲.ABAQUS動力無限元人工邊界研究[J].巖土力學(xué),2014,35(10):3 007-3 012.

[5]夏凡,董良國,馬在田.三維彈性波數(shù)值模擬中的吸收邊界條件[J].地球物理學(xué)報,2004,47(1):132-136.

[6]Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics.1985,50(9):705-708.

[7]Kosloff D, Kosloff R. Absorbing boundaries for wave propagation problems[J].Computer.Phys., 1986,63(6):363-376.

[8]Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics.1994,114(7):185-200.

[9]Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bulletin of the Seismological Society of America.1977,67(10):1 529-1 540.

[10]Higdon R L. Absorbing boundary conditions for elastic waves[J]. Geophysics.1991,56(2):231-241.

[11]Randall C J. Absorbing boundary condition for the elastic wave equation[J].Geophysics, 1988,53(8):611-624.

[12]Reynolds. Boundary conditions for the numerical solution of wave propagation problems[J]. Geophysics,1978,43(4):1099.

[13]R F Ungless. An infinite finite element[D]. Canada: University of British Columbia, 1973.

[14]P Bettess, More on infinite elements[J]. Int.J.Numer,Meth.Eng.1980,15(11):1 613-1 626.

[15]G Beer, J L M eek. Infinite domain element[J] Int.J.Numer,Meth.Eng.1981(17):43-52.

[16]O.C.Zienkiewicz et al. A novel boundary infinite element[J].Int.J.Numer,Meth.Eng.1983(19):393-404.

[17]賈坤,蔣樹屏,李建軍.無限單元在隧道工程數(shù)值分析中的應(yīng)用[J].西部探礦工程,2007(7):135-138.

[18]朱以文，蔡元奇，徐晗.ABAQUS與巖土工程分析[M].香港：中國圖書出版社，2005:136-137.

[19]費康，張建偉.ABAQUS在巖土工程中的應(yīng)用[M].北京：中國水利水電出版社，2013:54.

[20]夏媛媛,趙民,藏歌,等.正演模擬在解釋反演中的應(yīng)用[J].工程地球物理學(xué)報,2014,11(6):842-846.

[21]丘斌煌,晏紅艷,李三福,等.底辟模糊區(qū)地震資料處理關(guān)鍵技術(shù)[J].工程地球物理學(xué)報,2014,11(6):824-831.

[22]王敏.地震資料處理中的迭代法衰減多次波技術(shù)[J].工程地球物理學(xué)報,2015,12(6):800-807.

[23]楊凱,劉偉.基于希爾伯特—黃變換的地震資料高分辨率處理[J].工程地球物理學(xué)報,2015,12(1):22-26.

The Application of Infinite Element Boundaries to the Numerical Simulation of Finite Element Earthquake

Li Zhi1，Xue Tingxiao1,2，Zhang Hailiang3

(1.FacultyofEarthSciences,GuilinUniversityofTechnology,GuilinGuangxi541004,China; 2.KeyLaboratoryofGeologicalEngineeringCenterofGuangxiProvince,GuilinUniversityofTechnology,GuilinGuangxi541004,China； 3.ShandongGeologicalProspectingInstitute,ChinaChemicalGeologyandMineBureau,JinanShandong250013,China)

When the finite element method is used to research the numerical simulation of the earthquake, certain boundary conditions are need to apllied to the boundary of the model as for the problems in researching infinite zone. Other methods are more difficult to set up boundary conditions, but using the infinite element method can be more effective to design the absorbing boundary and programming is also more convenient. In this paper, some simple models are designed, the finite element method is used to simulate the model, and the infinite element is designed in the boundary of the model.Then the earthquake snapshot and the artificial earthquake record are obtained.The results of the analysis clearly show that the influence of the boundary of the model is almost eliminated after using the infinite element method, and it can also be used in complex terrain. In solving the problem of infinite domain, infinite element can truly simulate seismic wave propagation process.

seismic numerical modeling; finite element; infinite element

1672—7940(2016)01—0035—06

10.3969/j.issn.1672-7940.2016.01.006

桂林理工大學(xué)科研啟動項目(編號：002401003302)

李 志(1989-)，男，碩士研究生，主要從事地震數(shù)值模擬研究。E-mail：lentureli@163.com

薛霆虓(1974-)，男，副教授，主要從事地震勘探、地球動力學(xué)科研和教學(xué)工作。E-mail：xuetx@glite.edu.cn

P631.4

A

2015-07-06