李義豐，汪勝祥

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

基于計算機有限存儲空間和節(jié)約計算時間的考慮，在模擬計算聲波在無限或半無限介質(zhì)中的傳播時要解決由計算邊界的截斷而引起的邊界反射問題。因此，必須要有合適的邊界條件來減小或消除由于邊界截斷而引起的反射。完全匹配層Perfectly Matched Layer(PML)是由Berenger［1］于1994年提出的一種應(yīng)用于電磁波截斷計算的邊界吸收條件。其是一種魯棒且有效的吸收邊界條件［2］，理論證實在進行數(shù)值離散前，該PML層對任何頻率及任何方向上的入射波在截斷面上均不發(fā)生反射［1，3］。其作為經(jīng)典技術(shù)被廣泛地應(yīng)用于聲波［4］及彈性波［3，5］的數(shù)值計算中。但其缺點是要對計算場進行分裂，中間變量較多，場方程較為繁瑣，對傾斜入射波及表面波吸收效果較差。

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

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

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

在直角坐標系中，基于牛頓運動定律，二維吸收流體媒質(zhì)的時域聲波方程如下:

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

1.1 CPML頻域聲波方程

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

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

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

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

其中，

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

1.2 CPML時域聲波方程

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

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

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

其中:

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

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

2 數(shù)值仿真計算

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

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

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

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

運算時，kmax=0，αmax=2π ×105，n1=3，n2=0，n3=1和:

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

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

計算時，

上式(32)中，右邊第一項為式(12)，f'為所用信號源函數(shù)式(25)在頻域中的表達式。計算時頻率步長為Δf=20 kHz，頻率變化范圍為fmin=20 kHz，fmax=1 MHz。圖2為頻域中的計算結(jié)果，分別為波形在頻率為f=40 kHz，380 kHz，740 kHz和 960 kHz 時的波形快照圖。從中我們可以定性看出CPML匹配層對各個頻率的聲波起到了很好的吸收效果，在截斷邊界處沒有反射發(fā)生。

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

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

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

圖4 損耗流體及無損耗流體介質(zhì)的聲壓在A、B兩點的比較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兩點的8層CPML的相對誤差Fig.5 Relative error at points A and B for 8 cell CPML

圖6 A和B兩測試點6層CPML及10層CPML的相對誤差Fig.6 Relative error at points A and B for 6 cell and 10 cell CPMLs

圖7 A和B兩測試點10層CPML及PML的相對誤差Fig.7 Relative error at points A and B for 10 cell CPML and PML

3 結(jié)論

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

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

［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］李義豐，李國峰，王 云.卷積完全匹配層在兩維聲波有限元計算中的應(yīng)用［J］.聲學(xué)學(xué)報，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.