完全匹配層中衰減函數(shù)的參數(shù)優(yōu)化分析

周鳳璽,張家齊,張海威

(1.蘭州理工大學(xué)土木工程學(xué)院,甘肅蘭州730050;2.西部土木工程防災(zāi)減災(zāi)教育部工程研究中心,甘肅蘭州730050)

?

完全匹配層中衰減函數(shù)的參數(shù)優(yōu)化分析

周鳳璽1,2,張家齊1,張海威1

(1.蘭州理工大學(xué)土木工程學(xué)院,甘肅蘭州730050;2.西部土木工程防災(zāi)減災(zāi)教育部工程研究中心,甘肅蘭州730050)

完全匹配層(perfectly matched layer,PML)是一種高效的吸收邊界條件,對體波和面波都有非常好的吸收效果,被廣泛應(yīng)用于彈性波數(shù)值模擬。針對二維彈性動力學(xué)問題,基于PML最大反射系數(shù)的理論推導(dǎo),對PML衰減函數(shù)中的參數(shù)取值進(jìn)行了優(yōu)化分析。首先基于復(fù)伸展坐標(biāo)變換,給出了一種適用于二階彈性波動方程的非分裂吸收邊界條件;然后通過求解平面P-SV波的波數(shù),得到了非分裂完全匹配層反射系數(shù)的解析表達(dá)式;最后采用COLLINO給出的衰減函數(shù)形式,令PML最大反射系數(shù)為最小,得到了衰減函數(shù)中PML的厚度、理論反射系數(shù)以及最大反射系數(shù)之間的相互關(guān)系。通過數(shù)值算例分析了PML最大反射系數(shù)的變化規(guī)律,為PML參數(shù)的選擇和優(yōu)化提供了理論依據(jù)。

復(fù)伸展坐標(biāo)變換;數(shù)值模擬;波數(shù);衰減函數(shù);完全匹配層

在地球物理學(xué)和地震學(xué)等研究領(lǐng)域,完全匹配層(PML)被廣泛應(yīng)用于數(shù)值模擬無邊界問題。在對電磁波、聲波、彈性波等進(jìn)行數(shù)值模擬的過程中,受計算機內(nèi)存和計算時間的限制,通常會將無限區(qū)域截斷,從而產(chǎn)生無意義的人為邊界反射。為了解決這一問題,有關(guān)專家和學(xué)者提出了大量的吸收邊界條件,其中由BERENGER[1-2]提出的完全匹配層在理論上可以吸收來自各個方向、各種頻率的波,不產(chǎn)生任何邊界反射。

完全匹配層最初是BERENGER于1994年針對電磁波問題的研究提出的一種有效吸收邊界條件[1]。同年,CHEW等[3]將復(fù)伸展坐標(biāo)應(yīng)用在PML吸收邊界條件中。之后又有許多專家學(xué)者將完全匹配層問題的研究擴展到了聲波、彈性波等數(shù)值模擬中,如HASTINGS等[4]將PML應(yīng)用到一階速度-應(yīng)力聲波方程中;COLLINO等[5],CHEW等[6]以及裴正林[7]將PML應(yīng)用到彈性波的數(shù)值模擬中;DIMITRI等[8]提出了適用于二階彈性波動方程分裂式的PML吸收邊界條件;BECACHE等[9]在理論上證明了各向異性彈性波動方程PML吸收邊界條件的穩(wěn)定性。在完全匹配層問題的研究分析中,產(chǎn)生了多種形式的PML,主要有分裂式PML(Splitting PML,SPML)和非分裂式PML(Non-splitting PML,NPML)兩種,兩者除在計算的復(fù)雜程度和計算量上有區(qū)別外,對人工邊界反射有相同的吸收效果,在完全匹配層問題分析中都得到了廣泛應(yīng)用。

在地震波正演和波動方程偏移數(shù)值計算中,邊界處反射波的能量嚴(yán)重影響了計算結(jié)果的精確性。為了減小完全匹配層的邊界反射,BERENGER[1],HASTINGS等[4],COLLINO等[5]分別針對一階分裂Maxwell差分方程、一階分裂式速度-應(yīng)力聲波方程和一階分裂式速度-應(yīng)力波動方程進(jìn)行了完全匹配層問題的優(yōu)化;VADIM[10]針對各向同性彈性波動方程提出了最優(yōu)有限差分完全匹配層;COLLINO等[11]針對曲線坐標(biāo)系下的Maxwell方程提出了最優(yōu)完全匹配層;SHIMADA等[12]對有限元離散在完全匹配層邊界上所產(chǎn)生的反射能量進(jìn)行了研究。對于一階分裂式和非分裂式完全匹配層的優(yōu)化工作還有很多[13-19],但關(guān)于二階非分裂完全匹配層的優(yōu)化研究很少。本文針對二階彈性波動方程的非分裂PML吸收邊界條件[20-21],通過對最大反射系數(shù)的理論推導(dǎo),得出了該PML的反射機理,對于非分裂完全匹配層問題的研究有一定的參考意義。

1 基本方程

考慮二維彈性半空間動力學(xué)問題,用矩陣形式表示的彈性波基本方程如下:

(1)

σ=Cε

(2)

ε=Lu

(3)

其中,

式中:σ表示應(yīng)力張量;ε表示應(yīng)變張量;u表示位移矢量;f表示外力矢量;C表示彈性張量;λ和μ為拉梅系數(shù);ρ表示彈性介質(zhì)的密度;(¨)表示對時間的二階偏導(dǎo);L表示線性微分算子。

2 二階非分裂PML吸收邊界條件

考慮彈性波在半空間無限平面區(qū)域Ω=(-∞,∞)×[0,∞)內(nèi)傳播,可將無限區(qū)域進(jìn)行人為截斷,計算區(qū)域限制為一個有限的矩形區(qū)域,所有彈性波的波源限制在矩形計算區(qū)域ΩRD=(-a1,a1)×(0,a3),(a1,a3>0)的內(nèi)部。在ΩRD區(qū)域內(nèi),所有彈性波以任意角度向外傳播,在ΩRD區(qū)域外側(cè),考慮波的傳播速度不變,并在計算區(qū)域ΩRD的外圍布置厚度為L的完全匹配層,用于吸收傳出區(qū)域ΩRD的彈性波,如圖1所示。

圖1 半空間無限區(qū)域的完全匹配層二維截斷

(4a)

(4b)

(5)

對(4)式進(jìn)行Fourier變換,即將波動方程由時域變換到頻域,彈性波動方程可以改寫為:

(6a)

(6b)

引入復(fù)伸展坐標(biāo)以及連續(xù)的衰減函數(shù)ζi:

(7)

復(fù)伸展坐標(biāo)的偏微分形式為:

(8)

考慮在吸收層外時衰減函數(shù)ζi=0,在吸收層內(nèi)時,衰減函數(shù)ζi為正。

利用公式(7)對(6)式進(jìn)行復(fù)伸展坐標(biāo)變換,可以得到:

(9a)

(9b)令

(10)

按照GROTE等[22]所給出的求解過程,利用復(fù)伸展坐標(biāo)的偏微分公式(8)對(9)式進(jìn)行替換,并在(9)式兩側(cè)同時乘以γ1γ3,可以得到:

(11a)

(11b)

根據(jù)(10)式,可以得到以下關(guān)系式:

-ω2γ1γ3=-ω2+iω(ζ1+ζ3)+ζ1ζ3

(12a)

(12b)

(12c)

將(12)式代入(11)式可以得到:

(13a)

(13b)

(14a)

(14b)

(14c)

(14d)

也可寫成如下形式:

(15a)

(15b)

(15c)

(15d)

將(14)式代入(13)式可以得到:

(16b)

將(16)式由頻域轉(zhuǎn)換到時域,可以得到彈性波動方程的非分裂PML吸收邊界條件的表達(dá)式:

(17a)

(17b)

同時,將(15)式轉(zhuǎn)化到時間域,可以得到輔助函數(shù)的控制方程:

(18a)

(18b)

(18c)

(18d)

在區(qū)域ΩRD的內(nèi)部,衰減函數(shù)ζi(i=1,3)和輔助矢量Φ會消失,因此公式(17)還原為公式(4)。在彈性波動方程非分裂PML吸收邊界條件的推導(dǎo)過程中,只需引入4個標(biāo)量輔助函數(shù)φ1,φ2,φ3,φ4,即可將完全匹配層引入到彈性波動方程中,而且在推導(dǎo)過程中不需要提高時間和空間的階數(shù),形式簡單,易于實現(xiàn)。

3 完全匹配層反射系數(shù)的求解

基于二階彈性波動方程的非分裂PML吸收邊界條件,考慮彈性波沿著x1方向和x3方向呈指數(shù)遞減,設(shè)縱波和橫波的勢函數(shù)分量形式為[4]:

u1=D1ei(ωt-k1x1-k3x3)

(19a)

u3=D3ei(ωt-k1x1-k3x3)

(19b)

(19c)

(19d)

(19e)

(19f)

將(19)式代入(17)式可以得到方程:

ζ3)D1-ρζ1ζ3D1

(20a)

ζ3)D3-ρζ1ζ3D3

(20b)

再將(19)式代入(18)式可以得到位移和勢函數(shù)幅值之間的轉(zhuǎn)換關(guān)系:

(21a)

(21b)

(21c)

(21d)

將(21)式代入(20)式可以得到方程:

ω2ρD1-iωρ(ζ1+ζ3)D1-ρζ1ζ3D1

(22a)

ω2ρD3-ρ(ζ1+ζ3)iωD3-ρζ1ζ3D3

(22b)

將(22)式合并,利用Snell定律k1/k3=sinθ/cosθ(θ為彈性波的入射角),整理可得:

(23)

式中:σ=λ+2μ。(23)式經(jīng)過一系列的計算并整理后可得到兩個控制方程:

(24a)

(24b)

其中,參數(shù)a,b,c的表達(dá)式為:

(25a)

(25b)

c=[ω2-iω(ζ1+ζ3)-ζ1ζ3]ρ

(25c)

將(25)式代入(24)式,并以P波的波數(shù)求解為例,可以得到:

(26)

S波的波數(shù)求解過程與P波相同。為了簡化(26)式,令二維函數(shù)f(ζ1,ζ3)為:

(27)

應(yīng)用一階二維McLaughlin展式對(27)式連續(xù)函數(shù)進(jìn)行求解,取頻率項為正可以得到:

f(ζ1,ζ3)=ω-ζ1isin2θ-ζ3icos2θ

(28)

因此,完全匹配層中P波的波數(shù)表達(dá)式為:

(29a)

(29b)

S波的波數(shù)表達(dá)式為:

(30a)

(30b)

將(30a)式、(30b)式代入(19a)式、(19b)式可以得到該吸收邊界條件的位移表達(dá)式:

(31a)

(31b)

考慮衰減函數(shù)趨于0,(31)式可以還原為以波數(shù)表示的彈性波在單相固體介質(zhì)中傳播的位移表達(dá)式。

(32b)

式中:RP和RS分別表示完全匹配層中P波和S波的反射系數(shù)。

為進(jìn)一步求解反射系數(shù),采用COLLINO等[11]提出的衰減函數(shù)形式:

(33)

式中:L為完全匹配層厚度;r為理論反射系數(shù);v為彈性波的速度。

將(33)式代入(32)式任意一項可以得到:

(34)

完全匹配層的優(yōu)化即是令PML最大反射系數(shù)為最小。通過對(34)式分析可知,當(dāng)θ=0,π/2時,最大反射系數(shù)Rmax表達(dá)式為:

Rmax=eLlogr

(35)

分析(35)式可知,當(dāng)完全匹配層厚度為0,或者理論反射系數(shù)為1時,完全匹配層處于全反射狀態(tài),此時最大反射系數(shù)的值Rmax=1。

4 數(shù)值分析

為了更好地說明完全匹配層最大反射系數(shù)與完全匹配層厚度以及理論反射系數(shù)之間的關(guān)系,也為了在完全匹配層數(shù)值模擬實驗分析時能夠確定衰減函數(shù)中參數(shù)的最優(yōu)值,對(35)式進(jìn)行了數(shù)值計算分析。

1) 在不同厚度、不同理論反射系數(shù)的情況下,分析完全匹配層最大反射系數(shù)的變化。取完全匹配層厚度L在0～2.5m變化,理論反射系數(shù)r取10-j,j=1,2,3,4,5,6,由(35)式可以得到最大反射系數(shù)Rmax的變化曲線,如圖2所示。當(dāng)理論反射系數(shù)r一定時,完全匹配層的厚度L越大,最大反射系數(shù)越小;當(dāng)完全匹配層的厚度L一定時,理論反射系數(shù)r取值越小,最大反射系數(shù)越小。即理論反射系數(shù)取值越小,完全匹配層厚度越大時,完全匹配層的吸收效果越好。

2) 在不同厚度、不同最大反射系數(shù)的情況下,分析理論反射系數(shù)對數(shù)的變化?？紤]完全匹配層厚度L在0～10m變化,最大反射系數(shù)Rmax取值為10-j,j=1,2,3,4,5,6,由(35)式得到理論反射系數(shù)的對數(shù)logr的變化曲線,如圖3所示。當(dāng)最大反射系數(shù)Rmax一定時,完全匹配層的厚度L越大,理論反射系數(shù)的對數(shù)越大;當(dāng)完全匹配層的厚度L一定時,最大反射系數(shù)Rmax越大,理論反射系數(shù)的對數(shù)越大。據(jù)此性質(zhì)可對COLLINO等[11]提出的衰減函數(shù)中參數(shù)的取值進(jìn)行優(yōu)化。

圖2 最大反射系數(shù)Rmax與完全匹配層厚度L的關(guān)系曲線

圖3 理論反射系數(shù)的對數(shù)logr與完全匹配層厚度L的關(guān)系曲線

5 結(jié)束語

本文基于復(fù)伸展坐標(biāo)變換給出了一種適用于二階彈性波動方程的非分裂吸收邊界條件,理論上推導(dǎo)了該完全匹配層的最大反射系數(shù)解析表達(dá)式。通過對最大反射系數(shù)進(jìn)行理論分析和數(shù)值分析,得出以下結(jié)論:①最大反射系數(shù)隨著完全匹配層厚度的增加而減小,隨理論反射系數(shù)的減小而減小;②理論反射系數(shù)的對數(shù)隨完全匹配層厚度以及最大反射系數(shù)的增加而增加。研究結(jié)果對于完全匹配層數(shù)值模擬實驗分析中衰減函數(shù)的參數(shù)取值有一定的指導(dǎo)意義,當(dāng)衰減函數(shù)的參數(shù)取值適當(dāng)時,能夠消除完全匹配層的邊界反射。

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

[2] BERENGER J P.Three-dimensional perfectly mathced layers for the absorption of eletromagnetic waves[J].Journal of Computational Physics,1996,127(2):363-379

[3] CHEW W C,WEEDON W H.A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates[J].Microwave and Optical Technlogy Letters,1994,7(13):599-604

[4] HASTINGS F,SCHNEIDER J B,BROSCHAT S L.Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[J].Journal of the Acoustic Society of America,1996,100(5):3061-3069

[5] COLLINO F,TSOGKA C.Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heteregeneous media[J].Geophysics,2001,66(1):294-307

[6] CHEW W C,LIU Q H.Perfectly matched layers for elastodynamic:a new absorbing boundary condition[J].Journal of Computational Acoustics,1996,4(4):341-359

[7] 裴正林.三維各向同性介質(zhì)彈性波方程交錯網(wǎng)格高階有限差分法模擬[J].石油物探,2005,44(4):308-315 PEI Z L.Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method[J].Geophysical Prospecting for Petroleum,2005,44(4):308-315

[8] DIMITRI K,JEROEN T.A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J].Geophysical Journal International,2003,154(1):146-153

[9] BECACHE E,FAUQUEUX S,JOLY P.Stability of perfectly matched layers,group velocities and anisotropic waves[J].Journal of Computational Physics,2003,188(2):399-433

[10] VADIM L.Optimal discretization of PML for elasticity problems[J].Electronic Transactions on Numerical Analysis Etna,2008,30(7):258-277

[11] COLLINO F,MONK P B.Optimizing the perfectly matched layer[J].Computer Methods in Applied Mechanics and Engineering,1998,164(1):157-171

[12] SHIMADA T,HASEGAWA K,SATO S.Analysis of reflection powers from a perfectly matched layer for elastic waves in the frequency domain finite element model[J].Japanese Journal of Applied Physics,2011,50(7):913-919

[13] JOLY P.An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[J].SeMA Journal,2012,57(1):5-48

[14] 陳可洋.完全匹配層吸收邊界條件研究[J].石油物探,2010,49(5):472-477 CHEN K Y.Study on perfectly matched layer absorbing condition[J].Geophysical Prospecting for Petroleum,2010,49(5):472-477

[15] 陳可洋.聲波完全匹配層吸收邊界條件的改進(jìn)算法[J].石油物探,2009,48(1):76-79 CHEN K Y.Improved algorithm for absorbing boundary condition of acoustic perfectly matched layer[J].Geophysical Prospecting for Petroleum,2009,48(1):76-79

[16] 高剛,賀振華,黃德濟,等.完全匹配層人工邊界條件中的衰減函數(shù)分析[J].石油物探,2011,50(5):430-433 GAO G,HE Z H,HUANG D J,et al.Analysis on attenuation factor in the processing of artifical boundary conditions of PML[J].Geophysical Prospecting for Petroleum,2011,50(5):430-433

[17] 王永剛,邢文軍,謝萬學(xué),等.完全匹配層吸收邊界條

件的研究[J].中國石油大學(xué)學(xué)報(自然科學(xué)版),2007,31(1):19-24 WANG Y G,XING W J,XIE W X,et al.Study of absorbing condition by perfectly matched layer[J].Journal of China University of Petroleum (Edition of Natural Science),2007,31(1):19-24

[18] 熊章強,毛承英.聲波數(shù)值模擬中改進(jìn)的非分裂式PML邊界條件[J].石油地球物理勘探,2011,46(1):35-39 XIONG Z Q,MAO C Y.Improved unsplit PML boundary conditions in acoustic wave numerical simulation[J].Oil Geophysical Prospecting,2011,46(1):35-39

[19] 秦臻,任培罡,姚姚.彈性波正演模擬中PML吸收邊界條件的改進(jìn)[J].地球科學(xué)(中國地質(zhì)大學(xué)學(xué)報),2009,34(4):658-664 QIN Z,REN P G,YAO Y.Improvement of PML absorbing boundary conditions in elastic wave forward modeling[J].Earth Science (Journal of China University of Geosciences),2009,34(4):658-664

[20] 周鳳璽,曹小林,賈克M B.極坐標(biāo)系下非分裂PML及時域有限元實現(xiàn)[J].應(yīng)用數(shù)學(xué)和力學(xué),2015,36(9):956-969 ZHOU F X,CAO X L,JAKSA M B.A non-splitting PML for elastic waves in polar coordinates and its finite element implementation[J].Applied Mathematics and Mechanics,2015,36(9):956-969

[21] 周鳳璽,高貝貝.多孔介質(zhì)瞬態(tài)分析中非分裂PML及時域有限元實現(xiàn)[J].應(yīng)用數(shù)學(xué)和力學(xué),2016,37(2):195-209 ZHOU F X,GAO B B.A non-splitting PML for transient analysis of poroelastic media and its finite element implementation[J].Applied Mathematics and Mechanics,2016,37(2):195-209

[22] GROTE M J,SIM I.Efficient PML for the wave equation[EB/OL].[2016-01-20]http:∥arxiv.org/pdf/1001.0319v1.pdf

(編輯:戴春秋)

The parameter optimization analysis for the attenuation function in the perfectly matched layer

ZHOU Fengxi1,2,ZHANG Jiaqi1,ZHANG Haiwei1

(1.DepartmentofGeotechnicalEngineering,LanzhouUniversityofTechnology,Lanzhou730050,China;2.TheWesternCivilEngineeringDisasterPreventionandMitigationEngineeringResearchCenteroftheMinistryofEducation,Lanzhou730050,China)

Perfectly matched layer (PML) is a high-efficiency absorbing boundary condition (ABC) to body waves and surface waves,and is widely applied to the numerical simulation of the elastic wave.Considering the two-dimensional elastic dynamic problems,on the basis of the theoretical solution of PML’s maximum reflection coefficient,we carry out optimization analysis on the parameters selection for the attenuation function in the PML.Primarily,based on the complex-stretching-coordinate transform,we deduced one kind of the unsplit absorbing boundary condition appropriate for the second-order elastic wave equation.Then by solving the wave number of plane P-SV waves,we obtained the analytical expression of maximum reflection coefficient of the PML.We chose the attenuation function form proposed by Francis Collino,by ordering the maximum reflection coefficient approaching the minimum value,and eventually obtained the mutual relations of the PML’ s thickness,the theoretical reflection coefficient and the maximum reflection coefficient.Finally,we analyzed the variation rule of PML’ s maximum reflection coefficient by numerical testing,which provided theoretical guidance for the selection and optimization of PML’ s parameters.

complex-stretching-coordinate transform,numerical simulation,wave number,attenuation function,perfectly matched layer

2016-01-25;改回日期：2016-04-05。

周鳳璽(1979—),男,博士,教授,主要從事巖土力學(xué)和非均勻材料結(jié)構(gòu)力學(xué)方面的研究和教學(xué)工作。

國家自然科學(xué)基金(51368038)、甘肅省環(huán)保廳科研項目(GSEP-2014-23)和甘肅省教育廳研究生導(dǎo)師基金(1103-07)聯(lián)合資助。

This research is financially supported by the National Natural Science Foundation (Grant No.61368038),the Science Project of of Environmental Department in Gansu Province (Grant No.GSEP-2014-23) and the Tutor Funded Project of Education Department in Gansu Province (Grant No.1103-07).

P631

A

1000-1441(2016)06-0793-07

10.3969/j.issn.1000-1441.2016.06.003