基于CPML的FDTD法礦山地質(zhì)災(zāi)害應(yīng)急數(shù)值模擬

尤雙雙，謝杰文，彭谷香

（湖北省核工業(yè)地質(zhì)局，湖北 孝感 432000）

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

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

1 基本原理

1.1 CPML的背景

CPML的理論是基于坐標(biāo)伸縮Maxwell導(dǎo)出的平面波在分界面的無反射條件。

頻域中修正的無源Maxwell：

經(jīng)推導(dǎo)得伸縮媒質(zhì)中平面波的波阻抗為：

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

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

1.2 時(shí)域CPML的公式

即將含指數(shù)函數(shù)的卷積化為循環(huán)卷積，方便進(jìn)行FDTD步進(jìn)計(jì)算。

為了簡(jiǎn)單起見，對(duì)CPML中x方向的磁場(chǎng)表達(dá)式進(jìn)行FDTD差分離散，如下所示：

2 數(shù)值模擬

2.1 一維情況

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

圖1 時(shí)間步t=65時(shí)，電場(chǎng)的傳播情況

圖2 時(shí)間步t=65時(shí)，磁場(chǎng)的傳播情況

圖3 時(shí)間步t=100時(shí)，電場(chǎng)的傳播情況

圖4 時(shí)間步t=100時(shí)，磁場(chǎng)的傳播情況

圖5 時(shí)間步t=110時(shí)，電場(chǎng)的傳播情況

圖6 時(shí)間步t=110時(shí)，磁場(chǎng)的傳播情況

圖7 時(shí)間步t=140時(shí)，電場(chǎng)的傳播情況

圖8 時(shí)間步t=140時(shí)，磁場(chǎng)的傳播情況

2.2 二維情況

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

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

圖9 時(shí)間步t=240時(shí)的電場(chǎng)平均功率

圖10 時(shí)間步t=310時(shí)的電場(chǎng)平均功率

圖11 時(shí)間步t=400時(shí)的電場(chǎng)平均功率

圖12 時(shí)間步t=500時(shí)的電場(chǎng)平均功率

圖13 時(shí)間步t=120時(shí)，電場(chǎng)的傳播情況

圖14 時(shí)間步t=230時(shí)，電場(chǎng)的傳播情況

圖15 時(shí)間步t=350時(shí)，電場(chǎng)的傳播情況

圖16 時(shí)間步t=500時(shí)，電場(chǎng)的傳播情況

3 結(jié)論

一維情況中的圖1和圖2表明在時(shí)間步t=65時(shí)，電磁波仍停留在有效區(qū)域傳播，并未到達(dá)CPML吸收層；同時(shí)說明了該礦山地質(zhì)災(zāi)害應(yīng)急模型中設(shè)置的CPML吸收層不影響波在有效區(qū)域內(nèi)傳播。圖3至圖6呈現(xiàn)出電磁波在衰減吸收區(qū)域被CPML層迅速吸收，圖7和圖8表明電場(chǎng)和磁場(chǎng)幾乎被吸收完畢。說明CPML層的吸收效果很好，在一維模型中設(shè)置比較合理。

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