完全匹配层中衰减函数的参数优化分析

周凤玺,张家齐,张海威

(1.兰州理工大学土木工程学院,甘肃兰州730050;2.西部土木工程防灾减灾教育部工程研究中心,甘肃兰州730050)



完全匹配层中衰减函数的参数优化分析

周凤玺1,2,张家齐1,张海威1

(1.兰州理工大学土木工程学院,甘肃兰州730050;2.西部土木工程防灾减灾教育部工程研究中心,甘肃兰州730050)

完全匹配层(perfectly matched layer,PML)是一种高效的吸收边界条件,对体波和面波都有非常好的吸收效果,被广泛应用于弹性波数值模拟。针对二维弹性动力学问题,基于PML最大反射系数的理论推导,对PML衰减函数中的参数取值进行了优化分析。首先基于复伸展坐标变换,给出了一种适用于二阶弹性波动方程的非分裂吸收边界条件;然后通过求解平面P-SV波的波数,得到了非分裂完全匹配层反射系数的解析表达式;最后采用COLLINO给出的衰减函数形式,令PML最大反射系数为最小,得到了衰减函数中PML的厚度、理论反射系数以及最大反射系数之间的相互关系。通过数值算例分析了PML最大反射系数的变化规律,为PML参数的选择和优化提供了理论依据。

复伸展坐标变换;数值模拟;波数;衰减函数;完全匹配层

在地球物理学和地震学等研究领域,完全匹配层(PML)被广泛应用于数值模拟无边界问题。在对电磁波、声波、弹性波等进行数值模拟的过程中,受计算机内存和计算时间的限制,通常会将无限区域截断,从而产生无意义的人为边界反射。为了解决这一问题,有关专家和学者提出了大量的吸收边界条件,其中由BERENGER[1-2]提出的完全匹配层在理论上可以吸收来自各个方向、各种频率的波,不产生任何边界反射。

完全匹配层最初是BERENGER于1994年针对电磁波问题的研究提出的一种有效吸收边界条件[1]。同年,CHEW等[3]将复伸展坐标应用在PML吸收边界条件中。之后又有许多专家学者将完全匹配层问题的研究扩展到了声波、弹性波等数值模拟中,如HASTINGS等[4]将PML应用到一阶速度-应力声波方程中;COLLINO等[5],CHEW等[6]以及裴正林[7]将PML应用到弹性波的数值模拟中;DIMITRI等[8]提出了适用于二阶弹性波动方程分裂式的PML吸收边界条件;BECACHE等[9]在理论上证明了各向异性弹性波动方程PML吸收边界条件的稳定性。在完全匹配层问题的研究分析中,产生了多种形式的PML,主要有分裂式PML(Splitting PML,SPML)和非分裂式PML(Non-splitting PML,NPML)两种,两者除在计算的复杂程度和计算量上有区别外,对人工边界反射有相同的吸收效果,在完全匹配层问题分析中都得到了广泛应用。

在地震波正演和波动方程偏移数值计算中,边界处反射波的能量严重影响了计算结果的精确性。为了减小完全匹配层的边界反射,BERENGER[1],HASTINGS等[4],COLLINO等[5]分别针对一阶分裂Maxwell差分方程、一阶分裂式速度-应力声波方程和一阶分裂式速度-应力波动方程进行了完全匹配层问题的优化;VADIM[10]针对各向同性弹性波动方程提出了最优有限差分完全匹配层;COLLINO等[11]针对曲线坐标系下的Maxwell方程提出了最优完全匹配层;SHIMADA等[12]对有限元离散在完全匹配层边界上所产生的反射能量进行了研究。对于一阶分裂式和非分裂式完全匹配层的优化工作还有很多[13-19],但关于二阶非分裂完全匹配层的优化研究很少。本文针对二阶弹性波动方程的非分裂PML吸收边界条件[20-21],通过对最大反射系数的理论推导,得出了该PML的反射机理,对于非分裂完全匹配层问题的研究有一定的参考意义。

1 基本方程

考虑二维弹性半空间动力学问题,用矩阵形式表示的弹性波基本方程如下:

(1)

σ=Cε

(2)

ε=Lu

(3)

其中,

式中:σ表示应力张量;ε表示应变张量;u表示位移矢量;f表示外力矢量;C表示弹性张量;λ和μ为拉梅系数;ρ表示弹性介质的密度;(¨)表示对时间的二阶偏导;L表示线性微分算子。

2 二阶非分裂PML吸收边界条件

考虑弹性波在半空间无限平面区域Ω=(-∞,∞)×[0,∞)内传播,可将无限区域进行人为截断,计算区域限制为一个有限的矩形区域,所有弹性波的波源限制在矩形计算区域ΩRD=(-a1,a1)×(0,a3),(a1,a3>0)的内部。在ΩRD区域内,所有弹性波以任意角度向外传播,在ΩRD区域外侧,考虑波的传播速度不变,并在计算区域ΩRD的外围布置厚度为L的完全匹配层,用于吸收传出区域ΩRD的弹性波,如图1所示。

图1 半空间无限区域的完全匹配层二维截断

(4a)

(4b)

(5)

对(4)式进行Fourier变换,即将波动方程由时域变换到频域,弹性波动方程可以改写为:

(6a)

(6b)

引入复伸展坐标以及连续的衰减函数ζi:

(7)

复伸展坐标的偏微分形式为:

(8)

考虑在吸收层外时衰减函数ζi=0,在吸收层内时,衰减函数ζi为正。

利用公式(7)对(6)式进行复伸展坐标变换,可以得到:

(9a)

(9b)令

(10)

按照GROTE等[22]所给出的求解过程,利用复伸展坐标的偏微分公式(8)对(9)式进行替换,并在(9)式两侧同时乘以γ1γ3,可以得到:

(11a)

(11b)

根据(10)式,可以得到以下关系式:

-ω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)式由频域转换到时域,可以得到弹性波动方程的非分裂PML吸收边界条件的表达式:

(17a)

(17b)

同时,将(15)式转化到时间域,可以得到辅助函数的控制方程:

(18a)

(18b)

(18c)

(18d)

在区域ΩRD的内部,衰减函数ζi(i=1,3)和辅助矢量Φ会消失,因此公式(17)还原为公式(4)。在弹性波动方程非分裂PML吸收边界条件的推导过程中,只需引入4个标量辅助函数φ1,φ2,φ3,φ4,即可将完全匹配层引入到弹性波动方程中,而且在推导过程中不需要提高时间和空间的阶数,形式简单,易于实现。

3 完全匹配层反射系数的求解

基于二阶弹性波动方程的非分裂PML吸收边界条件,考虑弹性波沿着x1方向和x3方向呈指数递减,设纵波和横波的势函数分量形式为[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)式可以得到位移和势函数幅值之间的转换关系:

(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)式经过一系列的计算并整理后可得到两个控制方程:

(24a)

(24b)

其中,参数a,b,c的表达式为:

(25a)

(25b)

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

(25c)

将(25)式代入(24)式,并以P波的波数求解为例,可以得到:

(26)

S波的波数求解过程与P波相同。为了简化(26)式,令二维函数f(ζ1,ζ3)为:

(27)

应用一阶二维McLaughlin展式对(27)式连续函数进行求解,取频率项为正可以得到:

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

(28)

因此,完全匹配层中P波的波数表达式为:

(29a)

(29b)

S波的波数表达式为:

(30a)

(30b)

将(30a)式、(30b)式代入(19a)式、(19b)式可以得到该吸收边界条件的位移表达式:

(31a)

(31b)

考虑衰减函数趋于0,(31)式可以还原为以波数表示的弹性波在单相固体介质中传播的位移表达式。

(32b)

式中:RP和RS分别表示完全匹配层中P波和S波的反射系数。

为进一步求解反射系数,采用COLLINO等[11]提出的衰减函数形式:

(33)

式中:L为完全匹配层厚度;r为理论反射系数;v为弹性波的速度。

将(33)式代入(32)式任意一项可以得到:

(34)

完全匹配层的优化即是令PML最大反射系数为最小。通过对(34)式分析可知,当θ=0,π/2时,最大反射系数Rmax表达式为:

Rmax=eLlogr

(35)

分析(35)式可知,当完全匹配层厚度为0,或者理论反射系数为1时,完全匹配层处于全反射状态,此时最大反射系数的值Rmax=1。

4 数值分析

为了更好地说明完全匹配层最大反射系数与完全匹配层厚度以及理论反射系数之间的关系,也为了在完全匹配层数值模拟实验分析时能够确定衰减函数中参数的最优值,对(35)式进行了数值计算分析。

1) 在不同厚度、不同理论反射系数的情况下,分析完全匹配层最大反射系数的变化。取完全匹配层厚度L在0～2.5m变化,理论反射系数r取10-j,j=1,2,3,4,5,6,由(35)式可以得到最大反射系数Rmax的变化曲线,如图2所示。当理论反射系数r一定时,完全匹配层的厚度L越大,最大反射系数越小;当完全匹配层的厚度L一定时,理论反射系数r取值越小,最大反射系数越小。即理论反射系数取值越小,完全匹配层厚度越大时,完全匹配层的吸收效果越好。

2) 在不同厚度、不同最大反射系数的情况下,分析理论反射系数对数的变化。考虑完全匹配层厚度L在0～10m变化,最大反射系数Rmax取值为10-j,j=1,2,3,4,5,6,由(35)式得到理论反射系数的对数logr的变化曲线,如图3所示。当最大反射系数Rmax一定时,完全匹配层的厚度L越大,理论反射系数的对数越大;当完全匹配层的厚度L一定时,最大反射系数Rmax越大,理论反射系数的对数越大。据此性质可对COLLINO等[11]提出的衰减函数中参数的取值进行优化。

图2 最大反射系数Rmax与完全匹配层厚度L的关系曲线

图3 理论反射系数的对数logr与完全匹配层厚度L的关系曲线

5 结束语

本文基于复伸展坐标变换给出了一种适用于二阶弹性波动方程的非分裂吸收边界条件,理论上推导了该完全匹配层的最大反射系数解析表达式。通过对最大反射系数进行理论分析和数值分析,得出以下结论:①最大反射系数随着完全匹配层厚度的增加而减小,随理论反射系数的减小而减小;②理论反射系数的对数随完全匹配层厚度以及最大反射系数的增加而增加。研究结果对于完全匹配层数值模拟实验分析中衰减函数的参数取值有一定的指导意义,当衰减函数的参数取值适当时,能够消除完全匹配层的边界反射。

[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] 裴正林.三维各向同性介质弹性波方程交错网格高阶有限差分法模拟[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].石油物探,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] 高刚,贺振华,黄德济,等.完全匹配层人工边界条件中的衰减函数分析[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] 王永刚,邢文军,谢万学,等.完全匹配层吸收边界条

件的研究[J].中国石油大学学报(自然科学版),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] 熊章强,毛承英.声波数值模拟中改进的非分裂式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].地球科学(中国地质大学学报),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.极坐标系下非分裂PML及时域有限元实现[J].应用数学和力学,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] 周凤玺,高贝贝.多孔介质瞬态分析中非分裂PML及时域有限元实现[J].应用数学和力学,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—),男,博士,教授,主要从事岩土力学和非均匀材料结构力学方面的研究和教学工作。

国家自然科学基金(51368038)、甘肃省环保厅科研项目(GSEP-2014-23)和甘肃省教育厅研究生导师基金(1103-07)联合资助。

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