完全匹配層在時域有限元彈性波數值模擬中的應用

覃發兵,高志偉,解皓楠,徐振旺

(1.長江大學地球物理與石油資源學院,湖北武漢430100;2.長江大學管理學院,湖北荊州434023;3.中國石油天然氣股份有限公司新疆油田分公司采油一廠,新疆克拉瑪依834000;4.中國石油天然氣股份有限公司遼河油田分公司勘探開發研究院,遼寧盤錦124010)

高精度彈性波數值模擬是彈性波動方程成像和反演的核心步驟之一。彈性波數值模擬時,從計算可行性考慮,無界的地下介質需加人工邊界截斷為有界區域,但在人工邊界處會出現虛假反射,從而嚴重影響數值模擬的精度。因此,需要在人工邊界處進行特殊處理以消除虛假反射對數值模擬結果的影響。

根據地震波場數值模擬時是否在求解區域最外層附加額外的幾何空間,可將近四十年發展起來的人工邊界處理技術分為兩類:吸收邊界方法和吸收衰減層方法。吸收邊界方法是在計算區域最外層網格節點上引入額外變量壓制虛假反射。由于沒有增加計算區域的幾何空間,吸收邊界方法產生的額外計算量較少。經典的吸收邊界方法為CLAYTON等[1]于1977年提出的旁軸近似法,該方法能有效吸收垂直入射的波場能量,但對于大角度入射能量吸收效果較差,難以滿足精度要求。針對該問題,HIGDON[2]將旁軸近似法推廣為高階形式,使其能吸收更大范圍角度入射的波場能量,但高階形式的旁軸近似實現過程較復雜,且引入的額外計算量較大。不同于吸收邊界方法,吸收衰減層方法是在實際物理求解區域外額外設置多個衰減層以逐漸吸收傳入到衰減層中的波場能量。該方法顯然會引入額外的計算量,但往往能吸收較大范圍入射角的波場能量。1994年,BERENGER[3]在求解Maxwell方程組時提出了目前最經典的吸收衰減層方法——完全匹配層(PML)。該方法理論上能使任意角度任意頻率成分的邊界入射波的反射系數為零,并因此而得名。1996年,CHEW等[4]首次將PML推廣到彈性波方程的求解。早期PML需要分裂波場來求解,這對于大規模數值模擬是不容忽視的計算量。2007年,KOMATITSCH和MARTIN在地震波場模擬中發展了不需要分裂波場的卷積完全匹配層(CPML)[5],顯著提高了PML的計算效率及對近掠射的邊界入射波的吸收效果。隨后,CPML技術被迅速推廣到聲波方程模擬及更為復雜介質的數值模擬中[6-9]。近期,高正輝等[10]在2.5維聲波近似方程數值模擬中提出了一種改進的PML方法,該方法可使邊界入射波在PML吸收層中進行兩次能量衰減,從而進一步提高了PML的吸收效果。

然而,關于PML的應用研究大多在有限差分法數值模擬方面[10-15],這是因為PML在一階形式的波動方程求解中較容易加載。有限元法數值模擬主要基于二階波動方程,對于PML的應用存在一定的障礙。2003年,KOMATITSCH等[16]提出了一組適用于二階波動方程的PML方程組,并將其成功應用于譜元法地震波場數值模擬。盡管如此,關于PML在有限元法數值模擬中的應用研究仍然較少[17-18]。

本文基于LIU等[19]提出的PML在譜元法中的加載方法,推導了分裂式PML在二維彈性波時間域有限元數值模擬中的加載過程,說明了理論反射系數選取對PML吸收效果的影響,討論了狄利克雷邊界在PML外邊界層和自由邊界處的正確使用方式及其對PML吸收效果的影響。

1 PML的基本原理

頻率域各向異性彈性波動方程可表示為:

(1)

圖1 目標區域和PML邊界層區域圖示

(2)

依據鏈式法有:

(3)

?x表示對x坐標的偏導(下文同)。則:

(4)

其中

(5)

PML從純數學上講,是將對實坐標的偏導做如(4)式所示對復數坐標偏導的復數變換。對應地,PML的物理意義可借助平面波解來很好地理解。將(2)式代入平面波解:

(6)

2 含PML的有限元控制方程

2.1 含PML的彈性波方程

以二維各向同性完全彈性介質為例,其本構方程為:

(7)

(8)

將(7)式代入(8)式,有:

(9)

(10)

又知:

(11)

將(11)式代入(10)式得:

將位移場x分量Ux分裂為Ux=Ux1+Ux2+Ux3,并將其代入(12)式中的第一個方程,整理后有:

(13)

同理,將Uz分裂為Uz=Uz1+Uz2+Uz3,并將其代入(12)式中第二個方程,整理后有:

(14)

(13)式和(14)式中均含有高次項,故引入以下4個中間變量:

(15)

將(15)式代入(13)式和(14)式中,有:

(16)

和

(17)

對(16)式和(17)式進行時間傅里葉反變換,可得包含了PML的彈性波動方程:

(18)

和

(19)

式中:uj(j=x,z)表示對應Uj(j=x,z)的時間域位移場量,pxx、pxz、pzx、pzz為時間域中間變量。

2.2 彈性波方程的等效積分弱形式

設Ω為全區域,Γ為PML外邊界,nx為單位外法向量的x分量,nz為單位外法向量的z分量。采用Galerkin法,取試探函數φx和φz分別為位移x分量和z分量的變分,將(18)式和(19)式分別乘以φx和φz并求其加權積分。根據分部積分性質有:

(20)

和

φz+C12?xuxnzφz)dΓ-

dzuz2φzdΩ

(21)

(20)式和(21)式即為包含PML的彈性波方程等效積分弱形式。

2.3 空間離散及有限元控制方程

有限元的主要思想是將全區域Ω的求解轉變為一系列子區域Ωe(或單元區域)上的局部求解,將這些局部解組合即得到全局解。這里將全區域分解為若干線性三角形單元,即每個單元內有3個節點,進行雙線性插值,則單元內部的位移場可近似為:

(22)

(23)

其中

(24)

(25)

其中

(26)

符號“?”表示兩個向量之間對應元素相乘。綜上所述,可將(20)式和(21)式寫成有限元控制方程形式:

(27)

和

(28)

式中:M為全局質量矩陣;Kx1、Kx2、Kx3、Kx4、Kx5、Kz1、Kz2、Kz3、Kz4和Kz5均為全局剛度矩陣;Cx、Cz、Cxz、Cxx和Czz為全局阻尼矩陣。

3 有限元控制方程的時間離散

觀察(27)式和(28)式可見,加載PML的彈性波有限元控制方程有如下一般形式:

采用Euler方法進行時間離散求解類似(30)式的方程,可得到四個中間變量pxx、pxz、pzx和pzz的時間遞推格式:

(31)

其中,Δt為時間步長。對于類似(29)式的方程,可用中心差分法得到:

(32)

和

(33a)

(33b)

(33c)

(34)

4 模型算例

4.1 全空間均勻介質模型

為了說明理論反射系數的設定對PML吸收效果的影響,進行了二維彈性波均勻全空間數值模擬。圖2a為均勻各向同性彈性介質模型,密度為2000kg/m3,縱波速度為1800m/s,橫波速度為1100m/s,泊松比約為0.202。震源為正z方向激發的方向力源,子波函數為雷克子波,主頻為5Hz,位于(2000m,1000m)點處。采用對稱結構化線性三角形網格(圖2b),先將全區域用邊長為10m的正方形單元剖分,然后將每個正方形單元劃分成4個等腰直角三角形。在全空間條件下(即無面波情形),最短波長為1100/5/2.5=88.0m,三角形最大邊長為10m,則一個波長內的采樣點數至少為8.8。此外,三角形單元的高最小為5m,線性三角形單元CFL(courant-friedrichs-lewy)數約為0.67,時間步長應滿足穩定性條件Δt≤0.67×5/1800=1.86ms,本算例中取時間步長為1ms。

根據前人研究[5],阻尼函數dx可表示為:

(35)

圖2 均勻各向同性彈性介質模型(a)與對稱結構化線性三角形網格(b)

4.2 半空間均勻介質模型

狄利克雷邊界條件不同于自由邊界條件,是一種剛性邊界條件,即波傳播至此會被完全反射回去。假設PML吸收層最外層區域為Ωu,那么狄利克雷邊界條件可表示為:

(36)

下面討論半空間均勻介質情況下,在PML吸收層最外層施加狄利克雷邊界條件對PML吸收效果的影響。震源位置設為(2000m,50m),z=0處為自由表面。圖5和圖6分別為未加載狄利克雷邊界條件和不包括地表自由表面部分的PML最外吸收層上加載了狄利克雷邊界條件的波場快照。從圖5a和圖6a 均可看到自由表面處有明顯的面波產生。對比圖5b和圖6b可見,PML最外吸收層上不加載狄利克雷邊界條件會降低PML邊界的穩定性,產生可見的虛假波場(圖5b黃色虛線圈內)。因此,在半空間均勻介質情形下,PML吸收層最外層加載狄利克雷邊界條件能增強PML的數值穩定性。

圖3 不同理論反射系數1.5s時刻的總速度場波場快照對比a 理論反射系數R=10-2; b 理論反射系數R=10-3; c 理論反射系數R=10-4; d 理論反射系數R=10-5

圖4 不同理論反射系數條件下(1700m,1300m)點處地震記錄a 速度場x分量; b 速度場z分量

圖5 未加狄利克雷邊界條件的總速度場波場快照a 2.0s時刻; b 2.5s時刻

需要指出的是,在PML吸收層最外層加載狄利克雷邊界條件存在一個嚴重的數值陷阱。圖7為在包括地表自由表面部分的所有PML最外吸收層上加載了狄利克雷邊界條件模擬的波場快照。可以看出,在自由表面左右兩端產生了非常嚴重的面波虛假反射,嚴重干擾了目標區域內的波場求解,說明處于地表自由表面的PML吸收層最外層部分不可加載狄利克雷邊界條件。

圖6 加載了狄利克雷邊界條件的總速度場波場快照a 2.0s時刻; b 2.5s時刻

圖7 在處于自由表面上的PML吸收層最外層上加載狄利克雷邊界條件的總速度場波場快照a 2.0s時刻; b 2.5s時刻

4.3 起伏地表非均勻介質模型

為了檢驗PML吸收層中包含狄利克雷邊界在復雜介質條件下的吸收效果,采用如圖8所示的起伏地表非均勻介質模型進行了有限元數值模擬測試。模型x方向3000m,最大深度為1000m,地表最大高程為210 m。上層介質密度為2100kg/m3,縱波速度2800m/s,橫波速度1500m/s;下層介質密度為2250kg/m3,縱波速度3400m/s,橫波速度1800m/s。震源函數為雷克子波,主頻20Hz,位于(1500m,100m)點處。模型全區域被剖分為686025個非結構化線性三角形單元,其最大邊長約為5.4m,高最小為3.1m。由于穩定性要求,時間采樣間隔取為0.5ms。PML吸收層厚度設置為170m,對應于一個最大主波長(約20個網格點)。理論反射系數取10-6,以保證復雜模型中PML吸收層能充分吸收不同方向不同波長的能量。

圖8 起伏地表非均勻介質模型

圖9為0.4s、1.0s和2.0s時刻的總速度場波場快照。其中,圖9a、圖9c和圖9e為PML吸收層最外層不加狄利克雷邊界時的模擬結果;圖9b、圖9d和圖9f為PML吸收層最外層加載了狄利克雷邊界條件時的模擬結果。從圖9可以看出,在0.4s和1.0s時刻,PML吸收層最外層加載和不加載狄利克雷邊界條件得到的波場信息一致;但在2.0s時刻,PML吸收層最外層不加狄利克雷邊界條件時,模擬結果出現了不穩定現象(見圖9e黃色虛線框部分),這進一步說明了在PML吸收層最外層加載狄利克雷邊界條件能增強PML的數值穩定性。

圖9 不同時刻總速度場波場快照(從上到下波場快照時刻依次為0.4s、1.0s和2.0s)a,c,e PML吸收層最外層未加狄利克雷邊界條件; b,d,f PML 吸收層最外層加載狄利克雷邊界條件

4.4 討論

在全空間均勻介質算例中,PML吸收層厚度一定時,理論反射系數越小,PML吸收效果越好。由于PML是一種衰減性質的邊界條件,這類加衰減層的邊界條件經驗上遵循衰減層越厚吸收效果越好的規律,因此推斷理論反射系數一定時,PML厚度越大,PML吸收效果越好。總之,理論反射系數越小,PML吸收層越厚,PML吸收效果越好。針對不同復雜程度的速度模型,理論反射系數和PML吸收層厚度的選擇可根據數值試驗確定。觀察分析PML吸收層最外層加載狄利克雷邊界條件對PML吸收效果的影響,我們認為其與狄利克雷邊界本身性質有關。狄利克雷邊界條件是一種剛性邊界條件,在有限區域內向外傳播的波場經過PML衰減層逐步衰減,傳播到PML最外層時遇到狄利克雷邊界,未被PML層完全吸收的波場會被完全反射回PML衰減層,進行第二次指數衰減。因此,理論分析認為,在PML最外層添加狄利克雷邊界條件可使傳播到人工邊界處的波場經歷兩倍厚度的PML吸收層的指數衰減,減小了時間遞推過程中計算誤差的傳遞,從而增強了PML的數值穩定性。但在自由表面PML吸收層區域加載狄利克雷邊界條件時,產生了嚴重的面波虛假反射,這是因為自由表面的Rayleigh面波沿自由表面傳播至PML吸收層時,遇到狄利克雷邊界產生了全反射,繼續傳播至PML外邊界自由表面角落時,如同產生了二次震源,最終造成強虛假反射。

5 結束語

本文推導和實現了PML吸收邊界條件在二階彈性波方程有限元數值解法中的變分形式,通過全空間均勻介質數值實驗,測試了理論反射系數對PML吸收效果的影響,發現當PML吸收層厚度約為半個最大主波長、理論反射系數設置小于(或等于)10-5時,PML吸收效果達到最優。

半空間均勻介質數值實驗結果表明,在PML吸收層最外層加載狄利克雷邊界條件可提高PML的數值穩定性,但是處于自由表面的PML吸收層最外層部分不可加載狄利克雷邊界條件,否則會產生嚴重的虛假反射。

起伏地表復雜介質數值模擬結果進一步表明,在PML吸收層最外層加載狄利克雷邊界條件對PML數值穩定性改善效果顯著,說明狄利克雷邊界條件的引入可作為提高PML數值穩定性的有效手段之一。