UWR-FDTD礦井電磁波數(shù)值分析方法

張高敏，劉 毅，彭 銘

(中國(guó)礦業(yè)大學(xué)(北京) 機(jī)電與信息工程學(xué)院，北京 100083)

煤礦井下5G、WIFI6、人員定位和車輛定位等無線系統(tǒng)[1]的設(shè)計(jì)和布置，需進(jìn)行礦井電磁波分析，以擴(kuò)大無線基站信號(hào)的覆蓋范圍，提高無線通信質(zhì)量，節(jié)約建設(shè)成本。國(guó)內(nèi)外學(xué)者用多種方法對(duì)礦井電磁波傳播特性展開研究[2-6]，但這些方法難以分析巷道中防爆電氣開關(guān)、機(jī)電設(shè)備、金屬管道、動(dòng)力和通信線纜等對(duì)礦井電磁波傳播的影響[7]。時(shí)域有限差分(Finite-Difference Time-Domain，F(xiàn)DTD)是計(jì)算電磁學(xué)中求解麥克斯韋方程組的一種基本方法，該方法自身充分考慮了空間內(nèi)各種介質(zhì)對(duì)電磁波的反射、衍射、透射的影響以及不同類型電磁波之間的相互作用，使仿真環(huán)境更接近實(shí)際目標(biāo)環(huán)境，被廣泛用于電磁問題的數(shù)值仿真[8-10]。

作為一種全波數(shù)值分析方法，F(xiàn)DTD需要存儲(chǔ)全部Yee網(wǎng)格上的電磁場(chǎng)分量及中間變量，用于時(shí)間步迭代，大量的網(wǎng)格以及矩陣運(yùn)算對(duì)計(jì)算機(jī)內(nèi)存和算力提出了較高的要求。為了減少內(nèi)存使用量，提升FDTD的計(jì)算效率，擴(kuò)大FDTD的應(yīng)用范圍，大量學(xué)者對(duì)此問題進(jìn)行了深入研究。KATO等[11]通過降低計(jì)算空間維度來減少存儲(chǔ)需求，在二維球坐標(biāo)系下用FDTD計(jì)算閃電回?fù)敉ǖ喇a(chǎn)生的電場(chǎng)波形。MIGUEL等[12]從共享內(nèi)存、分布式內(nèi)存以及矢量化不同CPU體系結(jié)構(gòu)等3個(gè)方面研究FDTD并行化計(jì)算，得出增加分布式節(jié)點(diǎn)，可提高計(jì)算性能等結(jié)論。MAI等用傳統(tǒng)FDTD計(jì)算粗網(wǎng)格區(qū)域的場(chǎng)值，用弱條件穩(wěn)定時(shí)域有限差分(Weakly Conditionally Stable -FDTD)計(jì)算細(xì)網(wǎng)格，粗網(wǎng)格和細(xì)網(wǎng)格區(qū)域交界處采用了空間插值方法來減少仿真時(shí)間和內(nèi)存需求[13]，但過多的插值運(yùn)算會(huì)降低計(jì)算精度。何欣波等[14]將Unconditionally Stable FDTD與傳統(tǒng)FDTD相結(jié)合克服了CFL(Courant-Friedrichs-Lewy)條件限制，由粗網(wǎng)格決定時(shí)間步長(zhǎng)，在細(xì)網(wǎng)格區(qū)域及相鄰的粗網(wǎng)格區(qū)域計(jì)算特征值，通過減少特征值分解矩陣大小提高計(jì)算效率。DUAN等在電磁波的傳播方向上設(shè)置一個(gè)可移動(dòng)窗口，可移動(dòng)窗口內(nèi)又設(shè)置了一個(gè)小窗口用來存儲(chǔ)大部分電磁場(chǎng)分量，當(dāng)小窗口到達(dá)移動(dòng)窗口邊界時(shí)把小窗口數(shù)據(jù)復(fù)制到新的移動(dòng)窗口開始新區(qū)域的計(jì)算[15]。AHDAB等[16]為降低計(jì)算資源，把全部的障礙物設(shè)在激勵(lì)源附近區(qū)域，用FDTD計(jì)算出該區(qū)域內(nèi)電磁波相互影響的全部結(jié)果，遠(yuǎn)場(chǎng)電磁波則用拋物線方程計(jì)算，BAKIRTZIS等[17]則提出了一種FDTD和射線追蹤相結(jié)合的電磁波傳播方法來降低FDTD的算力需求。

FDTD方法可以分析巷道中防爆電氣開關(guān)、機(jī)電設(shè)備、金屬管道、動(dòng)力和通信線纜等對(duì)礦井電磁波傳輸影響。但采用FDTD計(jì)算數(shù)百米三維立體巷道電磁波傳輸衰減，對(duì)計(jì)算設(shè)備的內(nèi)存、算力都提出了非常高要求，制約了FDTD在礦井電磁波衰減仿真的應(yīng)用。針對(duì)煤礦井下巷道斷面小，軸向長(zhǎng)的特點(diǎn)，筆者提出了單位窗口循環(huán)使用的時(shí)域有限差分方法(Unit Window Recycling，UWR-FDTD)。該方法沿巷道軸向把整個(gè)巷道均勻劃分成多個(gè)大小相同的虛擬窗口，僅為一個(gè)與虛擬窗口大小相同的單位窗口分配內(nèi)存空間，所有虛擬窗口根據(jù)自身編號(hào)循環(huán)使用單位窗口內(nèi)存進(jìn)行FDTD迭代，每一個(gè)虛擬窗口迭代完成后，存儲(chǔ)單位窗口電磁場(chǎng)，用于合成整個(gè)巷道內(nèi)的電磁場(chǎng)。UWR-FDTD方法不改變FDTD的迭代方程，可根據(jù)計(jì)算機(jī)硬件設(shè)定單位窗口大小，劃分窗口個(gè)數(shù)，F(xiàn)DTD迭代過程僅在一個(gè)單位窗口的內(nèi)存空間內(nèi)進(jìn)行，大幅減少了傳統(tǒng)FDTD方法對(duì)內(nèi)存的需求，提高了計(jì)算效率。

1 UWR-FDTD建模

1.1 FDTD原理

麥克斯韋方程組揭示了電和磁之間的轉(zhuǎn)換規(guī)律，預(yù)測(cè)了電磁波的存在。FDTD直接求解麥克斯韋方程組研究電和磁之間的關(guān)系，成為計(jì)算電磁學(xué)中的一個(gè)重要數(shù)值計(jì)算方法[18]。FDTD把計(jì)算空間劃分成多個(gè)Yee網(wǎng)格，電場(chǎng)可以設(shè)置在Yee網(wǎng)格棱的中心，磁場(chǎng)在Yee網(wǎng)格面的中心，這樣每個(gè)電場(chǎng)分量周圍都環(huán)繞4個(gè)磁場(chǎng)分量，每個(gè)磁場(chǎng)分量也環(huán)繞著4個(gè)電場(chǎng)分量，與自然界中電場(chǎng)與磁場(chǎng)的物理結(jié)構(gòu)保持一致。由于電場(chǎng)與磁場(chǎng)在時(shí)間和空間上都交叉分布，迭代過程也很容易用差分格式實(shí)現(xiàn)。用麥克斯韋方程組的兩個(gè)旋度方程，在時(shí)間和空間上差分可以推導(dǎo)出FDTD的迭代方程[19]。

(1)

式(1)中的E和H在直角坐標(biāo)系下沿x，y，z三個(gè)方向可以拆分成關(guān)于Ex，Ey，Ez和Hx，Hy，Hz的6個(gè)標(biāo)量方程。設(shè)定空間離散步長(zhǎng)和時(shí)間離散步長(zhǎng)，把上述6個(gè)標(biāo)量方程中關(guān)于x，y，z，t的偏導(dǎo)數(shù)用中心差分近似后，即可得到電磁場(chǎng)在x，y，z三個(gè)方向分量上FDTD顯示表達(dá)式。FDTD的計(jì)算空間包括待仿真的全部空間以及空間內(nèi)的激勵(lì)源和散射體，空間離散步長(zhǎng)一般小于待研究頻域內(nèi)最短波長(zhǎng)的1/10[16]。受計(jì)算機(jī)內(nèi)存容量限制，在模擬無限大區(qū)域內(nèi)的電磁特性時(shí)，還需要在相應(yīng)的截?cái)噙吔缣幗o出吸收邊界條件[21]。

1.2 UWR-FDTD模型

對(duì)于斷面為拱形或馬蹄形的礦井巷道可以用面積等效為一個(gè)矩形巷道，故在直角坐標(biāo)系下設(shè)等效矩形巷道斷面寬和高的大小分別為lx，ly，巷道軸向長(zhǎng)lz，3個(gè)方向Yee網(wǎng)格空間步長(zhǎng)均為δ。UWR-FDTD仿真時(shí)，用激勵(lì)源模擬發(fā)射天線，可以把電磁波信號(hào)通過更新方程耦合到空間網(wǎng)格內(nèi)。激勵(lì)源所在的巷道斷面為軸向起始參考位置，一般情況下激勵(lì)源位于巷道斷面中心，激勵(lì)源在巷道斷面上的不同位置可以研究發(fā)射天線位置對(duì)礦井電磁波傳播特性的影響。沿巷道軸向把巷道劃分成多個(gè)長(zhǎng)度大小均為lwz并從1順序編號(hào)的虛擬窗口，虛擬窗口的斷面大小與巷道的斷面大小相同。單位窗口中l(wèi)wz的大小可根據(jù)計(jì)算機(jī)硬件性能和巷道軸向長(zhǎng)度lz的大小來設(shè)定，一旦確定lwz后對(duì)計(jì)算結(jié)果向上取整，得單位窗口的數(shù)量nw：

(2)

仿真時(shí)，整個(gè)計(jì)算空間只需為電場(chǎng)和磁場(chǎng)分量及中間變量分配一個(gè)單位窗口的存儲(chǔ)空間用于FDTD迭代過程、與電磁波傳播lwz距離所用時(shí)間Tu相關(guān)的緩存空間。實(shí)際計(jì)算的巷道空間在3個(gè)方向上的大小分別為lx,ly和lz×nw，整個(gè)FDTD迭代過程在一個(gè)單位窗口空間內(nèi)完成。電磁波在一個(gè)單位窗口內(nèi)沿z方向的傳播時(shí)間Tu為

(3)

式中，Δt為Yee網(wǎng)格的時(shí)間步長(zhǎng)。

電磁波從激勵(lì)源經(jīng)過Tu時(shí)間到達(dá)單位窗口z方向的最后一個(gè)Yee網(wǎng)格，在此之后需要存儲(chǔ)最后一個(gè)網(wǎng)格的電場(chǎng)和磁場(chǎng)分量作為下一個(gè)虛擬窗口迭代的激勵(lì)源。

針對(duì)三維礦井巷道特點(diǎn)提出的UWR-FDTD方法同樣適用于一維和二維空間。使用UWR-FDTD方法對(duì)一維空間仿真時(shí)，Yee網(wǎng)格可比作一個(gè)點(diǎn)的擴(kuò)展，需要按照時(shí)間步存儲(chǔ)最后一個(gè)點(diǎn)的各電場(chǎng)和磁場(chǎng)分量作為下一窗口迭代的子激勵(lì)源。若使用UWR-FDTD方法對(duì)二維空間進(jìn)行仿真時(shí)，Yee網(wǎng)格可比作一條線的擴(kuò)展，需要按照時(shí)間步存儲(chǔ)最后一條線的各電場(chǎng)和磁場(chǎng)分量作為下一窗口迭代的子激勵(lì)源。同理使用UWR-FDTD對(duì)三維空間進(jìn)行仿真，Yee網(wǎng)格可作為一個(gè)面的擴(kuò)展，需要按照時(shí)間步存儲(chǔ)最后一個(gè)面的各電場(chǎng)和磁場(chǎng)分量作為下一窗口迭代的子激勵(lì)源。電磁波在空氣介質(zhì)中以近似光速?gòu)募?lì)源到達(dá)第一個(gè)窗口最后一個(gè)Yee網(wǎng)格的時(shí)間為Tu，并且激勵(lì)源以固定頻率周期變化，計(jì)算空間內(nèi)各位置的電磁場(chǎng)也將呈周期性變化，在此可緩存單位窗口內(nèi)最后一個(gè)網(wǎng)格Tu/Δt個(gè)時(shí)間步的電磁場(chǎng)作為下一個(gè)虛擬窗口的激勵(lì)源。每一個(gè)單位窗口運(yùn)行完最后一個(gè)時(shí)間步長(zhǎng)后都把該窗口的電場(chǎng)各分量轉(zhuǎn)存到合成電場(chǎng)中。UWR-FDTD的基本原理如圖1所示。

圖1 UWR-FDTD原理Fig.1 Principle of UWR-FDTD

為避免FDTD多個(gè)時(shí)間步迭代后的仿真結(jié)果發(fā)散，空間離散步長(zhǎng)和時(shí)間離散步長(zhǎng)上需要滿足CFL穩(wěn)定性條件[19-20]：

(4)

若網(wǎng)格3個(gè)方向上的空間步長(zhǎng)均為δ，則式(4)可簡(jiǎn)化為

(5)

S=cΔt/δ

(6)

1.3 邊界處理

在直角坐標(biāo)系中用x，y，z三個(gè)方向電場(chǎng)矢量和(下文用x方向表示計(jì)算空間的寬度方向，y方向表示高度方向，z方向表示長(zhǎng)度方向)表示Yee網(wǎng)格的合成電場(chǎng)，該電場(chǎng)選在Yee網(wǎng)格的左下角。在兩個(gè)虛擬窗口或傳播介質(zhì)與PML吸收邊界的交界面上，若使用左下角合成電場(chǎng)表示邊界處的電場(chǎng)，用最后一個(gè)Yee網(wǎng)格的電場(chǎng)表示傳播介質(zhì)和PML交界處的電場(chǎng)，則存在一個(gè)Yee網(wǎng)格空間步長(zhǎng)的計(jì)算誤差，如圖2中用Em(1,1,nm)表示第nm個(gè)網(wǎng)格的電場(chǎng)情況。因此，需要對(duì)邊界處的電場(chǎng)進(jìn)一步處理，減小該誤差。邊界電場(chǎng)處理方法如圖2所示。

圖2 邊界電場(chǎng)處理方法Fig.2 Treatment method of boundary electric fields

圖2中傳播介質(zhì)空間和PML空間沿x,y方向都包含nxny個(gè)Yee網(wǎng)格，傳播介質(zhì)空間沿z方向分別被劃分為nm個(gè)網(wǎng)格，PML包含na個(gè)網(wǎng)格。傳播介質(zhì)空間和PML空間第1個(gè)Yee網(wǎng)格的電場(chǎng)分別為Em(1,1,1)和Ea(1,1,1)，傳播介質(zhì)中最后一個(gè)Yee網(wǎng)格電場(chǎng)為Em(1,1,nm)?？梢越魄蠼釫m(1,1,nm+1)的值當(dāng)作傳播介質(zhì)和PML交界面處的電場(chǎng)，以及進(jìn)入PML層的初始電場(chǎng)Ea(1,1,1)，實(shí)現(xiàn)減少1個(gè)空間網(wǎng)格步長(zhǎng)的誤差。FDTD迭代過程可得在兩種介質(zhì)交界面上x和y方向上的電場(chǎng)Emx(1,1,nm+1)和Emy(1,1,nm+1)，由于Emz(1,1,nm+1)在PML空間無法用傳播介質(zhì)空間的電磁場(chǎng)分量直接計(jì)算，并且電磁波進(jìn)入到PML空間后電磁場(chǎng)強(qiáng)度的衰減速率按照指數(shù)增加，在此采用傳播介質(zhì)中Emz(1,1,nm-1)和Emz(1,1,nm)的線性插值表示Emz(1,1,nm+1)?？傻迷诮唤缑嫔涎刈鴺?biāo)軸3個(gè)方向上的電場(chǎng)為

(7)

2 數(shù)值計(jì)算與分析

本文實(shí)驗(yàn)仿真軟件為Matlab R2021a，操作系統(tǒng)為Windows 10 專業(yè)版，CPU型號(hào)為Intel(R) Core(TM) i7-7700K CPU @ 4.20 GHz，系統(tǒng)內(nèi)存容量32 GB。在此設(shè)計(jì)不同的實(shí)驗(yàn)場(chǎng)景從計(jì)算精度與計(jì)算效率2個(gè)方面對(duì)比UWR-FDTD和傳統(tǒng)FDTD的數(shù)值計(jì)算性能，并研究UWR-FDTD中窗口大小與仿真時(shí)間之間的關(guān)系。

2.1 內(nèi)存空間

煤礦巷道壁材料主要為混凝土，其電磁參數(shù)變化不大，可以等效為各向同性的一個(gè)狹長(zhǎng)的計(jì)算空間。設(shè)巷道斷面大小為nxny個(gè)Yee網(wǎng)格，軸向長(zhǎng)為nz個(gè)Yee網(wǎng)格，所需電磁場(chǎng)各分量及中間變量個(gè)數(shù)為ne，則傳統(tǒng)FDTD所需內(nèi)存空間κFDTD為

κFDTD=nxnynzne

(8)

(9)

單位窗口內(nèi)循環(huán)計(jì)算出的電場(chǎng)分量最后要合成整個(gè)巷道的電場(chǎng)，因此還需要分配一個(gè)額外內(nèi)存空間存儲(chǔ)整個(gè)巷道內(nèi)的電場(chǎng)κa，大小為

κa=nxnynz

(10)

由式(9)和(10)可得UWR-FDTD所需的全部?jī)?nèi)存空間κUWR為

(11)

式(8)給出了傳統(tǒng)FDTD所需的內(nèi)存空間，用式(11)除以式(8)可得UWR-FDTD所需的內(nèi)存空間與傳統(tǒng)FDTD所需內(nèi)存空間的比值，經(jīng)過化簡(jiǎn)后為

(12)

由式(12)可知ζ的大小與劃分窗口的個(gè)數(shù)以及電磁場(chǎng)及中間變量的個(gè)數(shù)有關(guān)。

下面在一維二維三維空間內(nèi)分別分析。根據(jù)FDTD更新方程，一維計(jì)算空間需1個(gè)電場(chǎng)和1個(gè)磁場(chǎng)變量，用指數(shù)差分BerengerPML吸收邊界條件至少需4個(gè)系數(shù)變量，故ne≥6。二維計(jì)算空間需2個(gè)電場(chǎng)、2個(gè)磁場(chǎng)變量、至少4個(gè)參與FDTD迭代的系數(shù)變量，故二維計(jì)算空間中ne≥8。對(duì)于三維計(jì)算空間，電磁場(chǎng)將分裂成12個(gè)電磁場(chǎng)變量，以及至少4個(gè)系數(shù)變量，故三維計(jì)算空間中ne≥16。不同維度空間所需的變量個(gè)數(shù)不同，劃分不同窗口數(shù)量時(shí)UWR-FDTD方法和傳統(tǒng)FDTD方法所需內(nèi)存的比值見表1。

表1 不同維度空間下ζ的大小

從表1可知，與傳統(tǒng)FDTD方法相比，用UWR-FDTD方法劃分的窗口數(shù)目越多，仿真時(shí)所需的內(nèi)存空間越??；在三維計(jì)算空間內(nèi)使用UWR-FDTD方法求解電磁問題時(shí)，節(jié)省的內(nèi)存空間遠(yuǎn)大于一維和二維計(jì)算空間。

2.2 計(jì)算精度

FDTD仿真過程需要用與時(shí)間相關(guān)的函數(shù)作為激勵(lì)源引入到計(jì)算區(qū)域，常用周期時(shí)諧函數(shù)和脈沖函數(shù)作為時(shí)諧源和脈沖源。時(shí)諧源帶寬窄，一次可以仿真一個(gè)頻點(diǎn)的電磁波，但可以直觀的看出電磁波的傳播特性。脈沖源具有寬頻帶特性，使用脈沖源進(jìn)行一次FDTD仿真，可以得到較寬頻帶上的電磁響應(yīng)特性[22]。在此采用正弦波作為時(shí)諧源，微分高斯脈沖作為脈沖源對(duì)UWR-FDTD方法的性能進(jìn)行分析，2種激勵(lì)源時(shí)間變化的函數(shù)方程為

(13)

對(duì)于式(13)中的微分高斯脈沖，當(dāng)t0=τ且fmaxτ=2(fmax為微分高斯脈沖中待分析的最高電磁波頻率)，可以用穩(wěn)定因子S和單位波長(zhǎng)包含的Yee網(wǎng)格數(shù)量dλ把Δt隱式表示，微分高斯脈沖可以離散為

(14)

在一維空間內(nèi)待求計(jì)算空間長(zhǎng)度為5 m，窗口長(zhǎng)度為2 m，根據(jù)UWR-FDTD方法整個(gè)計(jì)算空間將被劃分為3個(gè)窗口，實(shí)際計(jì)算空間長(zhǎng)度為6 m。Yee網(wǎng)格空間步長(zhǎng)為載波波長(zhǎng)的1/20[23]，CFL穩(wěn)定因子S為0.5，采用10層指數(shù)差分Berenger PML吸收邊界條件。以激勵(lì)源位置為起始點(diǎn)，用900 MHz的正弦波激勵(lì)源和微分高斯脈沖激勵(lì)，在無耗介質(zhì)中迭代720時(shí)間步后電場(chǎng)的時(shí)域波形如圖3所示。由于迭代720時(shí)間步電場(chǎng)剛傳播到整個(gè)計(jì)算空間的最后一個(gè)網(wǎng)格，正弦波激勵(lì)源最右側(cè)電場(chǎng)還沒有完全穩(wěn)定，使用UWR-FDTD方法劃分3個(gè)窗口合并后的電場(chǎng)與傳統(tǒng)的FDTD方法計(jì)算出的電場(chǎng)一致性較好，在2，4 m等2個(gè)窗口的交界處合成電場(chǎng)較為平滑。同一計(jì)算空間微分高斯脈沖激勵(lì)源在無耗介質(zhì)中用UWR-FDTD和FDTD方法計(jì)算出的時(shí)域電場(chǎng)波形也幾乎完全一致。

圖3 一維無耗介質(zhì)中的時(shí)域電場(chǎng)Fig.3 Time domain electric field in one-dimensionallossless medium

同無耗介質(zhì)計(jì)算空間大小一樣，設(shè)有耗介質(zhì)電導(dǎo)率σ為0.36 S/m，900 MHz正弦波和微分高斯脈沖激勵(lì)源傳播6 m，劃分3個(gè)窗口的UWR-FDTD與FDTD方法計(jì)算出的時(shí)域電場(chǎng)波形也一致，如圖4所示。

圖4 一維有耗介質(zhì)中的時(shí)域電場(chǎng)Fig.4 Time domain electric field in one-dimensional loss medium

在長(zhǎng)10 m、寬4 m的二維無耗介質(zhì)計(jì)算空間中，設(shè)置UWR-FDTD單位窗口長(zhǎng)3 m，寬4 m，可知需要把計(jì)算空間劃分4個(gè)虛擬窗口，實(shí)際計(jì)算空間長(zhǎng)12 m、寬4 m。設(shè)正弦激勵(lì)源頻率為900 MHz，Yee網(wǎng)格空間步長(zhǎng)16.67 mm，CFL穩(wěn)定因子為0.5，時(shí)間步長(zhǎng)27.8 ps，窗口四周采用10層指數(shù)差分BerengerPML吸收邊界條件。激勵(lì)源位于計(jì)算空間寬度方向上的中心位置，從激勵(lì)源開始在整個(gè)二維無耗計(jì)算空間內(nèi)沿計(jì)算空間長(zhǎng)度方向的時(shí)域電場(chǎng)分布如圖5所示。圖5(a)，(b)為2種方法求解出在整個(gè)二維空間內(nèi)的瞬時(shí)電場(chǎng)強(qiáng)度分布情況，圖5(c)為接收點(diǎn)位于寬度方向上的中心位置，沿巷道長(zhǎng)度方向上各空間采樣點(diǎn)處的瞬時(shí)電場(chǎng)強(qiáng)度。

圖5 二維無耗介質(zhì)中的時(shí)域電場(chǎng)Fig.5 Time domain electric field in two-dimensionallossless medium

由圖5可以看出，UWR-FDTD方法劃分4個(gè)窗口單獨(dú)運(yùn)行后再合成的電場(chǎng)分布，與在一個(gè)窗口內(nèi)使用FDTD運(yùn)行的電場(chǎng)分布，在幅值和周期兩個(gè)方面都能較好吻合，表明UWR-FDTD能夠在二維計(jì)算空間內(nèi)取得較好計(jì)算結(jié)果。

用文獻(xiàn)[2]中高3.4 m、寬4.8 m的一個(gè)等效矩形巷道來仿真三維無耗計(jì)算空間，窗口四周采用10層指數(shù)差分BerengerPML吸收邊界條件。當(dāng)正弦激勵(lì)源頻率為900 MHz時(shí)，波長(zhǎng)為333.11 mm，CFL穩(wěn)定因子為0.40時(shí)，時(shí)間步長(zhǎng)為22.22 ps，Yee網(wǎng)格空間步長(zhǎng)為電磁波波長(zhǎng)的1/20，對(duì)應(yīng)單位窗口寬高分別288和204個(gè)Yee網(wǎng)格空間步。把偶極子激勵(lì)源放在巷道斷面中心，設(shè)單位窗口長(zhǎng)為300個(gè)Yee網(wǎng)格，巷道長(zhǎng)度約為5 m，用傳統(tǒng)FDTD與UWR-FDTD兩種方法，對(duì)3個(gè)虛擬窗口組成的巷道和5個(gè)虛擬窗口組成的巷道進(jìn)行仿真，接收位置同樣位于巷道斷面中心沿巷道軸向采樣，最終得到各采樣點(diǎn)處的合成電場(chǎng)時(shí)域波形如圖6所示。

由圖6可知，在三維巷道內(nèi)，用UWR-FDTD與傳統(tǒng)FDTD方法計(jì)算出的電場(chǎng)大小和周期較為一致，劃分多個(gè)虛擬窗口分別計(jì)算，最后合成的波形在2個(gè)窗口過渡時(shí)較為平滑，無明顯窗口分割痕跡。經(jīng)過在不同維度空間對(duì)UWR-FDTD方法計(jì)算精度的分析可知，該方法可用于煤礦井下、金屬和非金屬等其他礦井、鐵路隧道、公路隧道、地鐵、地下人防工程、城市地下管廊、地下商場(chǎng)等空間內(nèi)的電磁波分析。

圖6 三維無耗介質(zhì)中時(shí)域電場(chǎng)Fig.6 Time domain electric field in three-dimensionallossless medium

2.3 計(jì)算效率

UWR-FDTD方法把計(jì)算空間分割成多個(gè)虛擬窗口，電磁場(chǎng)分量及中間變量矩陣大小和單位窗口大小相同。傳統(tǒng)FDTD在1個(gè)時(shí)間步內(nèi)迭代整個(gè)計(jì)算空間的電磁場(chǎng)變量矩陣，兩者相比，UWR-FDTD在每個(gè)時(shí)間步的計(jì)算量大幅減小，理論上可以縮短仿真時(shí)間，提高計(jì)算效率，下面在三維巷道內(nèi)用仿真實(shí)驗(yàn)給予驗(yàn)證。

《煤礦安全規(guī)程》第九十條和《煤礦巷道斷面和交岔點(diǎn)設(shè)計(jì)規(guī)范》4.1.1節(jié)規(guī)定了軌道機(jī)車運(yùn)輸?shù)南锏纼舾咦攒壍赖烂嫫鸩坏玫陀? m，采(盤)區(qū)內(nèi)上山、下山和平巷的凈高度不得低于2m，薄煤層內(nèi)不得低于1.8 m，巷道凈寬不宜小于2.0 m。用900 MHz正弦波作為激勵(lì)源，穩(wěn)定因子S取0.4，空間步長(zhǎng)為電磁波波長(zhǎng)的1/20，約為16.67 mm。設(shè)巷道斷面寬高分別為140，110個(gè)Yee網(wǎng)格，約為2.33 m和1.84 m，單位窗口長(zhǎng)為140個(gè)Yee網(wǎng)格約為2.33 m。測(cè)試從1到22個(gè)單位窗口組成的三維巷道，長(zhǎng)度從2.33 m到51.33 m，在上述不同長(zhǎng)度的巷道內(nèi)用UWR-FDTD和FDTD兩種方法仿真同一電磁波傳播問題，仿真過程所需時(shí)間如圖7所示。

圖7 不同三維空間的仿真時(shí)間Fig.7 Simulation time in different three-dimensional spaces

當(dāng)單位窗口長(zhǎng)度和待仿真巷道長(zhǎng)度相等時(shí)，相當(dāng)于沒有劃分窗口，2種方法的仿真時(shí)間都為74 s，時(shí)間差為0，說明UWR-FDTD在不劃分窗口時(shí)與傳統(tǒng)FDTD的計(jì)算方法本質(zhì)相同。當(dāng)計(jì)算空間大小為8個(gè)虛擬窗口時(shí)，UWR-FDTD方法仿真用時(shí)2 739 s，F(xiàn)DTD方法仿真用時(shí)4 149 s，兩者用時(shí)相差1 410 s。當(dāng)計(jì)算空間長(zhǎng)度為22個(gè)虛擬窗口時(shí)，UWR-FDTD方法仿真用時(shí)19 499 s，F(xiàn)DTD方法仿真用時(shí)34 335 s，UWR-FDTD方法所需仿真時(shí)間比FDTD減少了14 836 s，降低了43.2%。故，UWR-FDTD方法非常適合狹長(zhǎng)空間的電磁問題數(shù)值仿真，計(jì)算空間越大，同F(xiàn)DTD相比UWR-FDTD方法的計(jì)算效率越高。

2.4 窗口大小優(yōu)化

對(duì)于同一段巷道，單位窗口的大小決定了需要?jiǎng)澐痔摂M窗口的數(shù)量，虛擬窗口數(shù)量又決定著每個(gè)虛擬窗口在單位窗口內(nèi)的迭代次數(shù)，迭代次數(shù)進(jìn)一步影響整個(gè)巷道的仿真時(shí)間。因此需要研究巷道斷面的大小、巷道長(zhǎng)度、電磁波頻率等多個(gè)條件下，虛擬窗口數(shù)量與仿真時(shí)間之間的關(guān)系，給出最優(yōu)單位窗口大小的計(jì)算方法。

2.4.1 巷道斷面

用900 MHz正弦激勵(lì)源模擬發(fā)射天線，激勵(lì)源位于巷道斷面中心，Yee網(wǎng)格空間步長(zhǎng)為波長(zhǎng)的1/20，巷道長(zhǎng)(z向)1 200個(gè)Yee網(wǎng)格，沿長(zhǎng)度方向分別劃分為2，3，4，5，6，8，10，12，15，20，24個(gè)窗口，用xy表示巷道寬和高的Yee網(wǎng)格數(shù)量。當(dāng)xy分別為140×110,90×70和50×40時(shí)，仿真上述11種虛擬窗口劃分方法所需的仿真時(shí)間見表2。

表2 不同巷道斷面大小的仿真時(shí)間

取表2中每一種xy仿真時(shí)間最小的一個(gè)虛擬窗口為時(shí)間基準(zhǔn)(加粗?jǐn)?shù)據(jù))，其他虛擬窗口的仿真時(shí)間都減去基準(zhǔn)虛擬窗口的仿真時(shí)間，可得虛擬窗口數(shù)量、xy大小與仿真時(shí)間之間的關(guān)系，如圖8所示。

由圖8可以看出，對(duì)于同一長(zhǎng)度的巷道，不管xy大小如何，劃分虛擬窗口數(shù)量過多或過少都會(huì)增加仿真時(shí)間，并且xy越大，劃分虛擬窗口的數(shù)量對(duì)仿真時(shí)間的影響越大。長(zhǎng)度為1 200個(gè)Yee網(wǎng)格的巷道劃分8～10虛擬窗口時(shí)所用的仿真時(shí)間相對(duì)較少，此時(shí)單位窗口的長(zhǎng)度為總巷道長(zhǎng)度的1/10～1/8。

2.4.2 巷道長(zhǎng)度

激勵(lì)源、穩(wěn)定因子和吸收邊界條件與2.4.1節(jié)中設(shè)定的數(shù)值相同，巷道斷面固定為140×110個(gè)Yee網(wǎng)格，巷道長(zhǎng)度分別為1 200，840和720個(gè)Yee網(wǎng)格。每種長(zhǎng)度巷道依舊劃分為2，3，4，5，6，8，10，12，15，20，24個(gè)虛擬窗口，使用UWR-FDTD方法仿真上述不同長(zhǎng)度巷道所需的仿真時(shí)間見表3。

表3 不同巷道長(zhǎng)度的仿真時(shí)間

由表3可知，巷道長(zhǎng)度方向Yee網(wǎng)格大小為1 200和840時(shí)，把上述巷道劃分為8個(gè)和5個(gè)虛擬窗口所需仿真時(shí)間最短；長(zhǎng)度為720個(gè)Yee網(wǎng)格，劃分6個(gè)虛擬窗口時(shí)所需仿真時(shí)間最短。根據(jù)表3數(shù)據(jù)，以仿真時(shí)間最短的虛擬窗口為基準(zhǔn)(加粗?jǐn)?shù)據(jù))，其他虛擬窗口的仿真時(shí)間都減去基準(zhǔn)窗口的仿真時(shí)間，可得不同巷道長(zhǎng)度與仿真時(shí)間的關(guān)系，如圖9所示。

圖9 巷道長(zhǎng)度與仿真時(shí)間的關(guān)系Fig.9 Relationship betweentunnel length and simulation time

由圖9可知不管待計(jì)算的巷道有多長(zhǎng)，劃分虛擬窗口數(shù)量過多或過少都會(huì)增加仿真時(shí)間。對(duì)于同一xy斷面，單位窗口長(zhǎng)度為對(duì)應(yīng)巷道長(zhǎng)度的1/10～1/6時(shí)，所需仿真用時(shí)相對(duì)較小。

2.4.3 電磁波頻率

用1 200，900，700 MHz的正弦波作為激勵(lì)源在長(zhǎng)×寬×高為840×90×70個(gè)Yee網(wǎng)格的巷道內(nèi)仿真，網(wǎng)格步長(zhǎng)都為電磁波波長(zhǎng)的1/10，計(jì)算不同電磁波頻率在劃分不同窗口數(shù)量時(shí)所需的仿真時(shí)間，統(tǒng)計(jì)劃分窗口數(shù)量與電磁波頻率之間的關(guān)系。3種電磁波頻率劃分不同窗口數(shù)量時(shí)所用的仿真時(shí)間見表4。

表4 窗口數(shù)量與不同頻率的仿真時(shí)間

每個(gè)頻率中選擇用時(shí)最少的虛擬窗口作基準(zhǔn)(加粗?jǐn)?shù)據(jù))，該頻率中其他數(shù)目虛擬窗口的仿真時(shí)間減去基準(zhǔn)窗口的時(shí)間，可得不同電磁波頻率對(duì)仿真時(shí)間的影響程度，如圖10所示。3條不同電磁波頻率的時(shí)間差曲線幾乎完全重合，說明UWR-FDTD方法在相同數(shù)量網(wǎng)格內(nèi)所需的仿真時(shí)間與電磁波頻率無關(guān)，但劃分虛擬窗口的數(shù)量與仿真時(shí)間關(guān)系較大。當(dāng)計(jì)算空間被分為8個(gè)虛擬窗口時(shí)，3種頻率都能在最短的時(shí)間內(nèi)完成數(shù)值仿真，此時(shí)單位窗口長(zhǎng)度為計(jì)算空間總長(zhǎng)度的1/8。

圖10 窗口數(shù)量和頻率與仿真時(shí)間的關(guān)系Fig.10 Relationship between electromagnetic wavefrequencies and simulation time

綜上所述，從巷道空間分析最優(yōu)單位窗口的劃分方法，由圖8可知,不同巷道斷面單位窗口長(zhǎng)度為巷道總長(zhǎng)度的1/10～1/6內(nèi)時(shí)，仿真時(shí)間最短；由圖9可知對(duì)于不同巷道長(zhǎng)度，單位窗口長(zhǎng)度也是在上述取值范圍內(nèi)仿真時(shí)間最短?？傻贸鲈谕挥?jì)算空間內(nèi)，單位窗口長(zhǎng)度為計(jì)算空間總長(zhǎng)度的1/10～1/6時(shí)仿真計(jì)算用時(shí)最短。對(duì)于不同電磁波頻率，從圖10中可得當(dāng)單位窗口長(zhǎng)度為計(jì)算空間總長(zhǎng)度的1/10～1/6時(shí)，計(jì)算效率較高。因此，在同一巷道內(nèi)，不管電磁波頻率多大，當(dāng)單位窗口長(zhǎng)度為計(jì)算空間總長(zhǎng)度的1/10～1/6時(shí)，都能取得較好的計(jì)算效率。

3 結(jié) 論

(1)在一維二維和三維計(jì)算空間內(nèi)，不管是無耗介質(zhì)還是有耗介質(zhì)，UWR-FDTD計(jì)算出的數(shù)值結(jié)果與傳統(tǒng)FDTD計(jì)算出的數(shù)值結(jié)果都一致，具有較高的計(jì)算精度。

(2)同傳統(tǒng)FDTD方法相比，用UWR-FDTD方法劃分的窗口數(shù)目越多，仿真時(shí)所需的內(nèi)存空間越??；在三維巷道中UWR-FDTD方法節(jié)省的內(nèi)存空間遠(yuǎn)大于一維和二維空間。

(3)同一巷道內(nèi)UWR-FDTD所需的仿真時(shí)間與劃分虛擬窗口的數(shù)量、巷道斷面面積、巷道長(zhǎng)度、電磁波頻率有關(guān)。單位窗口長(zhǎng)度為待計(jì)算巷道總長(zhǎng)度的1/10～1/6時(shí)，UWR-FDTD方法所需的仿真時(shí)間最短。