基于CPML的FDTD法礦山地質災害應急數值模擬

尤雙雙，謝杰文，彭谷香

（湖北省核工業地質局，湖北 孝感 432000）

PML是由Berenger[1]于1994年針對電磁波提出的一種非物理吸收邊界條件。因其在截斷面上具有對任意方向及任何頻率的波都不產生反射[2]的優勢而備受關注與運用。比如Chen等[3]在1997年采用PML吸收邊界條件求解電磁波方程，取得了不錯的成效；朱章虎等[4]學者于2006實現了PML-FDTD的編程；在礦震等方面也頗受好評，2007年，單啟銅等[5]研究了PML吸收邊界條件下的彈性介質波，取得了一系列的成果。但是需要分裂計算場使得其方程冗長而繁瑣[6]，后來Kuzuoglu[7]等在1996年提出了可以吸收低頻凋落模的復頻率參數完全匹配層(CFS-PML)。因為CFSPML吸收邊界條件不容易實現，所以在2000年時，Roden and Gedney[8]基于對坐標伸縮因子進行改進，從而提出了吸收效果好、易擴展等優點[9]的卷積完全匹配層(CPML)。葛德彪和閆玉波老師在著作[10]中談到，對于處理棱邊和角頂等區域，CPML比Gendey S.D.[11]（1996）提出的UPML更為簡潔和高效。由于CPML是由方程的連續性得到的，并沒有進行非物理分裂，因此其表達式的形式簡單，節省計算過程中的存儲量。在2005年，高春霞等[12]在模擬空中核爆電磁脈沖時推導了雙曲正交坐標系下的CPML截斷邊界的電磁場離散方程，王粵等[13]基于傅里葉變換的卷積定理，推導了三維金屬矩形波導在直角坐標系下的CPML的迭代公式。張顯文等[14]于2009年結合CFS-PML與高階差分法研究了彈性波方程。胡媛等[15]在2011年基于CUDA架構，采用并行算法實現了CPML-3D-FDTD。在2012年，朱翼超等[16]采用CPML吸收邊界條件模擬了脈沖輻射，而李義豐等[17]證明了CPML吸收邊界對聲波中的損耗體介質吸收果很好。李建等[18,19]于2014年驗證了CPML在2維TE波及3維HIE-FDTD的吸收性能很好。劉廣東[20,21]分別于2014年研究了Davidson-Cole介質的CPML以及2015年模擬了一般色散介質中的ADE-FDTD-CPML電磁波傳播情況。為了提供關于井下地質災害的理論依據，彭凌星等[22]在2017年利用CPML吸收邊界數值模擬了井下地質雷達探測，并且討論了關于復雜地質礦體的探測。

本文模擬了一維情況下的礦震脈沖波和二維情況下的TM波在CPML吸收邊界條件下的傳播情況。

1 基本原理

1.1 CPML的背景

CPML的理論是基于坐標伸縮Maxwell導出的平面波在分界面的無反射條件。

頻域中修正的無源Maxwell：

經推導得伸縮媒質中平面波的波阻抗為：

說明平面波的波阻抗不受伸縮因子的影響。

平面波在分界面處無反射條件為：表明若分界面兩側的橫向伸縮因子及媒質參數相等，則平面波在分界面處無反射。

1.2 時域CPML的公式

即將含指數函數的卷積化為循環卷積，方便進行FDTD步進計算。

為了簡單起見，對CPML中x方向的磁場表達式進行FDTD差分離散，如下所示：

2 數值模擬

2.1 一維情況

此模型設總網格數為200，主網格數為120，空間增量Δz=0.01m，時間步長為Δt=Δz c。帶寬為6的脈沖源位于計算域的中點，沿z方向在相對介電常數為1的均勻介質中傳播。兩端的吸收邊界均為CPML，其層數都為40，程序運行時間步數為140。以下各圖分別為各個時間步數時電場與磁場分量的傳播情況。

圖1 時間步t=65時，電場的傳播情況

圖2 時間步t=65時，磁場的傳播情況

圖3 時間步t=100時，電場的傳播情況

圖4 時間步t=100時，磁場的傳播情況

圖5 時間步t=110時，電場的傳播情況

圖6 時間步t=110時，磁場的傳播情況

圖7 時間步t=140時，電場的傳播情況

圖8 時間步t=140時，磁場的傳播情況

2.2 二維情況

該模型采用笛卡爾直角坐標系，自左至右為x軸，自下而上為y軸。

假設總網格數為340＊340，主網格數為300＊300，吸收邊界為CPML，層數為20。TM波的波源位于x軸中點的上方，源激發的波長λ為5.5× 10-7m，頻率為f=c0λ，相對介電常數為1。每個方向的單元格的大小為λ20。程序運行的時間步數為500，描繪了在各時間步時電場的平均功率及傳播情況。

圖9 時間步t=240時的電場平均功率

圖10 時間步t=310時的電場平均功率

圖11 時間步t=400時的電場平均功率

圖12 時間步t=500時的電場平均功率

圖13 時間步t=120時，電場的傳播情況

圖14 時間步t=230時，電場的傳播情況

圖15 時間步t=350時，電場的傳播情況

圖16 時間步t=500時，電場的傳播情況

3 結論

一維情況中的圖1和圖2表明在時間步t=65時，電磁波仍停留在有效區域傳播，并未到達CPML吸收層；同時說明了該礦山地質災害應急模型中設置的CPML吸收層不影響波在有效區域內傳播。圖3至圖6呈現出電磁波在衰減吸收區域被CPML層迅速吸收，圖7和圖8表明電場和磁場幾乎被吸收完畢。說明CPML層的吸收效果很好，在一維模型中設置比較合理。

對于二維TM波的傳播情況而言，我們應該結合電場的平均功率圖與電場的等值線圖進行相應的分析。由圖11和圖14、15可知，在時間步t=230之前，TM波仍然在均勻介質中傳播，還沒有到達吸收層。由圖11和圖12知，電場在CPML中衰減特別迅速；從圖14和圖15中規則的等值線可得，CPML的設置比較好，吸收反射波的效果不錯。