卷積完全匹配層在二維聲波吸收流體介質(zhì)中的有限元計(jì)算

李義豐，汪勝祥

(1.南京工業(yè)大學(xué) 電子與信息工程學(xué)院，南京 210009;2.法國(guó)電子、微電子及納米技術(shù)研究所，里爾 59652)

基于計(jì)算機(jī)有限存儲(chǔ)空間和節(jié)約計(jì)算時(shí)間的考慮，在模擬計(jì)算聲波在無限或半無限介質(zhì)中的傳播時(shí)要解決由計(jì)算邊界的截?cái)喽鸬倪吔绶瓷鋯栴}。因此，必須要有合適的邊界條件來減小或消除由于邊界截?cái)喽鸬姆瓷?。完全匹配層Perfectly Matched Layer(PML)是由Berenger［1］于1994年提出的一種應(yīng)用于電磁波截?cái)嘤?jì)算的邊界吸收條件。其是一種魯棒且有效的吸收邊界條件［2］，理論證實(shí)在進(jìn)行數(shù)值離散前，該P(yáng)ML層對(duì)任何頻率及任何方向上的入射波在截?cái)嗝嫔暇话l(fā)生反射［1，3］。其作為經(jīng)典技術(shù)被廣泛地應(yīng)用于聲波［4］及彈性波［3，5］的數(shù)值計(jì)算中。但其缺點(diǎn)是要對(duì)計(jì)算場(chǎng)進(jìn)行分裂，中間變量較多，場(chǎng)方程較為繁瑣，對(duì)傾斜入射波及表面波吸收效果較差。

現(xiàn)在，由Roden和 Gedney［6］所提出的基于復(fù)坐標(biāo)變換并利用復(fù)頻移擴(kuò)展坐標(biāo)形式吸收參數(shù)的卷積完全匹配層Convolution Perfectly Matched Layer(CPML)在數(shù)值計(jì)算中得到了廣泛的應(yīng)用，其被有效地應(yīng)用于不均勻介質(zhì)、有耗介質(zhì)、色散介質(zhì)、各向異性介質(zhì)及非線性介質(zhì)的截?cái)嘤?jì)算中［7-9］。與經(jīng)典PML相比，CPML的主要優(yōu)勢(shì)在于其不需要把場(chǎng)分裂開，形式簡(jiǎn)單，易于實(shí)現(xiàn);其對(duì)傾斜入射波和表面波有著比PML更好的吸收效果［10-12］;并有效地解決了漸消波邊界反射問題以及長(zhǎng)時(shí)間計(jì)算［13］而引起的邊界反射問題。

本文將CPML吸收邊界層應(yīng)用于二維吸收聲波的截?cái)嘤?jì)算中，分別在時(shí)域及頻域中推導(dǎo)出引入CPML邊界條件的吸收流體的波動(dòng)方程，并在有限元(FEM)計(jì)算軟件COMSOL中完成其時(shí)域的計(jì)算。計(jì)算過程中，通過在COMSOL中選擇與所推導(dǎo)出來的方程相對(duì)應(yīng)的計(jì)算模型，并設(shè)置相應(yīng)的系數(shù)參數(shù)來完成運(yùn)算。相對(duì)于CPML的時(shí)域有限差分 (FDTD)計(jì)算而言，在COMSOL中的有限元 (FEM)計(jì)算過程變得更加簡(jiǎn)單和易于實(shí)現(xiàn)。

1 吸收流體二維CPML聲波方程的推導(dǎo)

在直角坐標(biāo)系中，基于牛頓運(yùn)動(dòng)定律，二維吸收流體媒質(zhì)的時(shí)域聲波方程如下:

式中v、w分別為x、y方向上的粒子速度，p代表的是聲壓，ρ0、c0及r分別為媒質(zhì)的密度、聲速及吸收系數(shù)。吸收系數(shù)r用來描述物理介質(zhì)的吸收損失特性。通常對(duì)于吸收介質(zhì)而言，速度c及系數(shù)r是頻率依賴的。然而，在一定的頻率范圍內(nèi)，我們可以認(rèn)為c及r是可以頻率獨(dú)立的，在很多實(shí)驗(yàn)測(cè)量中這一點(diǎn)已經(jīng)被證實(shí)［14］。因此，為了簡(jiǎn)便起見，在本文的研究中介質(zhì)假定為是非頻散的，即系數(shù)c及r在所感興趣的頻帶范圍內(nèi)是頻率獨(dú)立的。

1.1 CPML頻域聲波方程

為了將CPML卷積完全匹配層引入到方程中，首先，將上述方程式(1)～(3)進(jìn)行傅里葉變換，得到頻域中的聲波方程。其次，在頻域中采用如下的復(fù)坐標(biāo)變換［15］:

式中，sx和sy是復(fù)頻移 (Complex Frequency Shifthed)擴(kuò)展坐標(biāo)變量，其表達(dá)式由Kuzuoglu和 Mittra［16］以如下形式給出:

其中，kx，y≥1，σx，y≥0 和 αx，y≥0。在 CPML 匹配層中，實(shí)數(shù)部分kx，y是為了衰減漸消波而引入的尺度因子，虛數(shù)部分 σx，y是該匹配層衰減透射波的衰減參數(shù)，αx，y是吸收層中完成Butterworth濾波的一個(gè)頻率依賴項(xiàng)［2］。在非 CPML 吸收層中，kx，y=1 且 σx，y=0。擴(kuò)展坐標(biāo)變量的使用使得CPML層具有比PML層更優(yōu)越的性能:它能更有效吸收表面波［12］，并有效解決由于大尺寸計(jì)算［11］及長(zhǎng)時(shí)間計(jì)算［13］所引起的計(jì)算不穩(wěn)定問題。這樣方程式(1)～(3)將變?yōu)槿缦滦问降姆匠探M:

式中為變量u的傅里葉變換。為了在有限元軟件COMSOL中實(shí)現(xiàn)對(duì)上述方程組的求解計(jì)算，方程式(8)～(10)需要被改寫成通用形式的穩(wěn)態(tài)偏微分(PDE)方程，最終得到了頻域形式的CPML方程為:

其中，

在sx和sy中，kx、σx及 αx的表達(dá)式將在數(shù)值仿真部分給出，ky、σy及 αy分別具有與kx、σx和 αx一樣的表達(dá)式，僅僅需要將x替換為y。

1.2 CPML時(shí)域聲波方程

為了得到時(shí)域中吸收聲波的CPML方程，需將頻域中CPML的方程式 (8)～(10)進(jìn)行傅里葉逆變換，并利用傅里葉變換的卷積定理，首先得到時(shí)域中卷積形式的聲波方程。卷積過程的消除過程可以參見文獻(xiàn)［17］，最終得到的消除了卷積且?guī)в杏洃涀兞浚?8-19］時(shí)間偏導(dǎo)的系統(tǒng)時(shí)域方程如下:

式中，δx= σx，/kx，βx= σx，/kx+ αx，δy和 βy分別是與 δx和βx類似的表達(dá)式，僅需用y替換x;A，B，C和D是為了消除卷積而引入的記憶變量，在非CPML吸收層，他們的值為零，所以他們的一階導(dǎo)數(shù)僅僅在整個(gè)計(jì)算區(qū)域的一小部分被求解。為了在COMSOL軟件中實(shí)現(xiàn)對(duì)上述方程組的求解，將方程式 (15)對(duì)時(shí)間再次求偏導(dǎo)后，把式(13)和式 (14)帶入其中便得到如下的兩階方程:

最后得到如下的在方程右邊僅含源項(xiàng)的系統(tǒng)方程:

其中:

式(22)中，k'x，y=kx，y-1，C'= ?C/?t，D'= ?D/?t。變量C'及D'對(duì)時(shí)間的偏導(dǎo)由如下兩個(gè)方程分別給出:

由式(21)可以看出，有關(guān)CPML的參數(shù)項(xiàng)只是作為源項(xiàng)f出現(xiàn)在最初的初始方程的右邊。相對(duì)于經(jīng)典PML層而言，CPML層中為了消除卷積而引入的記憶變量的數(shù)量(4個(gè))要少于PML層中由于需要把場(chǎng)分裂開而引入的變量的數(shù)量(6個(gè))，因而CPML層所需的總體存儲(chǔ)空間會(huì)變的更少。

2 數(shù)值仿真計(jì)算

在此部分，我們測(cè)試CPML匹配層對(duì)二維聲波在吸收流體介質(zhì)中的吸收性能，為此建立了如圖1所示的數(shù)值計(jì)算模型，19.2 mm×19.2 mm二維計(jì)算空間被劃分為64×64個(gè)FEM計(jì)算網(wǎng)格，在四周有8個(gè)計(jì)算網(wǎng)格厚度的CPML吸收層。信號(hào)源S位于數(shù)值計(jì)算空間的中心，選取A及B兩點(diǎn)作為聲波壓力場(chǎng)p的測(cè)試點(diǎn)。如圖所示，A點(diǎn)與B點(diǎn)距離左邊CPML層均為0.3 mm，A點(diǎn)位于與信號(hào)源S同一高度的水平線上，B點(diǎn)距離上邊的CPML層為0.3 mm。

圖1 CPML二維數(shù)值計(jì)算模型Fig.1 2D CPML numerical calculation model

數(shù)值計(jì)算所應(yīng)用的介質(zhì)參數(shù)為:聲速為c0=1 500 m/s，密度為 ρ0=1 000 kg/m3，吸收系數(shù)為r=2 ×10-5s/m2，最大的計(jì)算頻率為fmax=106Hz。所使用的信號(hào)源S是Gaussian函數(shù)的一階時(shí)間導(dǎo)數(shù)信號(hào)源，其具有1.5 mm的直徑，表達(dá)式為:

式中:tw=1.0 ×10-6s，信號(hào)源滯后時(shí)間t0=5.0 ×10-6s。在擴(kuò)展坐標(biāo)變量sx中，所應(yīng)用的x方向的CPML參數(shù)的表達(dá)式如下:

運(yùn)算時(shí)，kmax=0，αmax=2π ×105，n1=3，n2=0，n3=1和:

式中，x0和d分別為CPML層的起始位置坐標(biāo)和厚度;R0=2×10-5是理論上的反射系數(shù)。

首先，實(shí)現(xiàn)頻域的運(yùn)算，即在COMSOL軟件中利用其通用形式的PDE頻域計(jì)算模型式(30)對(duì)式(11)和式(12)進(jìn)行有限元求解。

計(jì)算時(shí)，

上式(32)中，右邊第一項(xiàng)為式(12)，f'為所用信號(hào)源函數(shù)式(25)在頻域中的表達(dá)式。計(jì)算時(shí)頻率步長(zhǎng)為Δf=20 kHz，頻率變化范圍為fmin=20 kHz，fmax=1 MHz。圖2為頻域中的計(jì)算結(jié)果，分別為波形在頻率為f=40 kHz，380 kHz，740 kHz和 960 kHz 時(shí)的波形快照?qǐng)D。從中我們可以定性看出CPML匹配層對(duì)各個(gè)頻率的聲波起到了很好的吸收效果，在截?cái)噙吔缣帥]有反射發(fā)生。

其次，進(jìn)行時(shí)域的求解，在COMSOL軟件中對(duì)時(shí)域CPML方程式(21)和式(22)進(jìn)行有限元求解。為了節(jié)約求解過程中的計(jì)算時(shí)間和存儲(chǔ)空間，將計(jì)算區(qū)域分為CPML區(qū)域和非CPML區(qū)域 (模擬區(qū)域);模擬區(qū)域使用常規(guī)的波動(dòng)方程，即f=0的方程式 (21);CPML區(qū)域使用CPML波動(dòng)方程，即帶有記憶變量A、B、C'、D'和源項(xiàng)f的方程式(21)和式(24)，在此A、B、C'、D'的時(shí)間偏導(dǎo)方程式(16)、式(17)、式(23)、式(24)一同被求解;將不同區(qū)域計(jì)算方程的各項(xiàng)分別與COMSOL軟件所提供的系數(shù)形式的PDE時(shí)域計(jì)算模型式(33)中的對(duì)應(yīng)項(xiàng)進(jìn)行對(duì)比，從而完成對(duì)模型中不同求解項(xiàng)系數(shù)矩陣(ea質(zhì)量系數(shù)、da衰減系數(shù)、c散射系數(shù)、α對(duì)流流量守恒系數(shù)、γ守恒流量源項(xiàng)、a吸收系數(shù)、β對(duì)流系數(shù))及f源項(xiàng)矩陣的設(shè)定，進(jìn)而完成對(duì)方程的有限元求解計(jì)算。圖3所示是針對(duì)于無損流體介質(zhì)的時(shí)域計(jì)算結(jié)果，為波形在不同時(shí)刻t=3.3 μs，6.3 μs，9.3 μs，12.3 μs，15.3 μs 和 27.3 μs 的快照?qǐng)D，從中我們可以定性的看出CPML起到了很好的吸收效果，在邊界上沒有發(fā)生聲波的反射。圖3(a)～(f)相對(duì)于圖3(b)的相對(duì)最大幅值為分別為 -1.48 dB，0 dB，-3.48 dB，-11.33 dB，-23.78 dB 和 -60.50 dB，從中我們可以看出聲波進(jìn)入CPML吸收層后，聲場(chǎng)能量被不斷的吸收。圖4為測(cè)試點(diǎn)A和B的聲場(chǎng)在流體介質(zhì)有損耗及無損耗情況下的聲壓比較，從中可以看出介質(zhì)的損耗使得A、B兩點(diǎn)的聲壓分別衰減9.08 dB和12.59 dB。

為了進(jìn)一步定量分析CPML對(duì)損耗流體介質(zhì)的吸收效果，在此計(jì)算了相對(duì)誤差，其定義為:

為了探討CPML層的厚度對(duì)波的吸收效果的影響，在此又計(jì)算了CPML分別為6層及10層時(shí)相對(duì)誤差隨時(shí)間變化的規(guī)律，如圖6所示。其中，當(dāng)CPML為6個(gè)網(wǎng)格厚度時(shí)，A、B兩點(diǎn)的最大相對(duì)誤差分別為-66.40 dB 及 -57.72 dB，當(dāng) CPML為10個(gè)網(wǎng)格厚度時(shí)，A、B點(diǎn)的最大相對(duì)誤差分別為-73.74 dB及-69.05 dB。從而可知增加CPML層的厚度，可以更好地吸收進(jìn)入其中的聲場(chǎng)能量，尤其是更好地吸收傾斜入射的聲波能量。為了比較CPML及PML層對(duì)入射波吸收效果的不同，在此對(duì)10個(gè)網(wǎng)格厚度的CPML及PML層的相對(duì)誤差進(jìn)行計(jì)算，所得A點(diǎn)的最大相對(duì)誤差分別為-73.74 dB和53.98 dB，B點(diǎn)的最大相對(duì)誤差分別為-69.05 dB和 -47.47 dB，誤差波形變化曲線如圖7所示。由此可以看出，無論對(duì)垂直入射波還是傾斜入射波，CPML有著比PML更好的吸收效果。

圖4 損耗流體及無損耗流體介質(zhì)的聲壓在A、B兩點(diǎn)的比較Fig 4 Results of comparison for the lossy and lossless media about the pressure field at the points A and B

圖5 損耗流體介質(zhì)中A及B兩點(diǎn)的8層CPML的相對(duì)誤差Fig.5 Relative error at points A and B for 8 cell CPML

圖6 A和B兩測(cè)試點(diǎn)6層CPML及10層CPML的相對(duì)誤差Fig.6 Relative error at points A and B for 6 cell and 10 cell CPMLs

圖7 A和B兩測(cè)試點(diǎn)10層CPML及PML的相對(duì)誤差Fig.7 Relative error at points A and B for 10 cell CPML and PML

3 結(jié)論

本文討論了卷積完全匹配層(CPML)在吸收流體介質(zhì)二維聲波方程數(shù)值計(jì)算中的應(yīng)用，推導(dǎo)給出了二維聲波CPML方程在頻域和時(shí)域的不同表達(dá)形式，并在有限元計(jì)算軟件COMSOL中完成對(duì)時(shí)域方程的數(shù)值求解，給出了取得較好吸收效果的各項(xiàng)參數(shù)的數(shù)值。

二維數(shù)值計(jì)算的結(jié)果表明對(duì)吸收流體介質(zhì)而言，卷積完全匹配層對(duì)進(jìn)入其中的聲場(chǎng)能量起到了很好的吸收效果。對(duì)8個(gè)網(wǎng)格厚度的CPML層來說，-71.23 dB及-66.10 dB的垂直入射和傾斜入射的相對(duì)誤差，說明其極大地衰減了進(jìn)入其中的能量，有效地消除了發(fā)生在邊界處的反射。通過計(jì)算比較6層CPML及10層CPML的相對(duì)誤差，可知增加CPML的厚度可以更加有效地衰減傾斜入射的聲場(chǎng)能量;同時(shí)，通過比較可知CPML比PML有著更好的吸收效果，能夠更有效地衰減進(jìn)入聲場(chǎng)的能量。

［1］Berenger J P.A perfectly matched layer for the absorption of electromagnetic waves［J］.J.Comput.Phys.，1994，114:195-200.

［2］Gedney S D.Perfeclty matched layer absorbing boundary conditions［C］//Computational Electrodynamics:The Finite-Difference Time-Domain Method.edited by A.Taflove and S.C.Hagness.Norwood:Artech House.2005，273-328.

［3］Chew W C， Liu Q H. Perfectlymatched layersfor elastodynamics:A new absorbing boundary condition［J］.J.Comp.Acoust.，1996，4(4):341-359.

［4］Liu Q H，Tao J.The perfectly matched layer for acoustic waves in absorptive media［J］.J.Acoust.Soc.Am.，1997，102(4):2072-2082.

［5］Hasting F D，Schneider J B，Broschat S L.Application of the perfectly matched layer(pml)absorbing boundary condition to elastic wave propagation［J］.J.Acoust.Soc.Am.，1996，100(5):3061-3069.

［6］Roden J A，Gedney S D.Convolution PML(CPML):An efficient FDTD implementation of the CFS-PML for arbitrary media［J］.Microwave Opt.Technol.Lett.，2000，27:334-339.

［7］Teixeira FL， Chew W C. Generalclosed-form PML constitutive tensorstomatch arbitrarybianisotropicand dispersive linear media［J］.IEEE Microwave and Guided Wave Letters.，1998，8:223-225.

［8］Bécache E，Petropoulos P G，Gedney S D.On the long-time behavior of unsplit perfectly matched layers［J］.IEEE Trans.Antennas Propagat.，2004，52(5):1335-1342.

［9］ Appel? D，Kreiss G.Application of the perfectly matched layer to the nonlinear wave equation［J］.Wave Motion，2007，44:531-548.

［10］Bou Matar O，Preobrazhensky V，Pernod P.Two-dimensional axisymetric numericalsimulation ofsupercriticalphase conjugation of ultrasound in active solid media［J］.J.Acoust.Soc.Am.，2005，118(5):2880-2890.

［11］ Drossaert F H，Giannopoulos A.Complex frequency shifted convolution PML for FDTD modeling of elastic waves［J］.Wave Motion，2007，44(7-8):593-604.

［12］ Komatitsch D，Martin R.An unsplit convolutional perfectly matched Layer improved at grazing incidence for the seismic wave equation［J］.Geophysics，2007，72(5):SM155-SM167.

［13］ Bécache E，Petropoulos P G，Gedney S D.On the long-time behavior of unsplit perfectly matched layers［J］.IEEE Trans.Antennas Propagat.，2004，52(5):1335-1342.

［14］Liu Q H，Chang C.Compressional head waves in attenuative formations:forward modelling and inversion ［J］.Geophysics.，1996，61(6):1908-1920.

［15］ Chew W C，Jin J M，Michelssen E.Complex coordinate stretching as a generalized absorbing boundary condition［J］.Microwave Opt.Technol.Lett.，1997，15:363-369.

［16］ Kuzuoglu M， MittraR. Frequencydependenceofthe constitutive parameters of causal perfectly matched anisotropic absorbers［J］.IEEE Microwave Guided Wave Lett.，1996，6:447-449.

［17］李義豐，李國(guó)峰，王 云.卷積完全匹配層在兩維聲波有限元計(jì)算中的應(yīng)用［J］.聲學(xué)學(xué)報(bào)，2010，35(6):601-607.

［18］Li Y F，BouMatar O.Convolutional perfectly matched layer for elastic second-order wave equation［J］.J.Acoust.Soc.Am.，2010，127(3):1318-1327.

［19］ Li Y F，Bou Matar O，Preobrazhensky V，et al.Convolution-Perfectly MatchedLayer(C-PML)absorbingboundary condition for wave propagation in piezoelectric solid［J］.In the Proceedings of 2008 IEEE Ultrasonic Symp.，2008，1568-1571.