張高敏,劉 毅,彭 銘
(中國礦業大學(北京) 機電與信息工程學院,北京 100083)
煤礦井下5G、WIFI6、人員定位和車輛定位等無線系統[1]的設計和布置,需進行礦井電磁波分析,以擴大無線基站信號的覆蓋范圍,提高無線通信質量,節約建設成本。國內外學者用多種方法對礦井電磁波傳播特性展開研究[2-6],但這些方法難以分析巷道中防爆電氣開關、機電設備、金屬管道、動力和通信線纜等對礦井電磁波傳播的影響[7]。時域有限差分(Finite-Difference Time-Domain,FDTD)是計算電磁學中求解麥克斯韋方程組的一種基本方法,該方法自身充分考慮了空間內各種介質對電磁波的反射、衍射、透射的影響以及不同類型電磁波之間的相互作用,使仿真環境更接近實際目標環境,被廣泛用于電磁問題的數值仿真[8-10]。
作為一種全波數值分析方法,FDTD需要存儲全部Yee網格上的電磁場分量及中間變量,用于時間步迭代,大量的網格以及矩陣運算對計算機內存和算力提出了較高的要求。為了減少內存使用量,提升FDTD的計算效率,擴大FDTD的應用范圍,大量學者對此問題進行了深入研究。KATO等[11]通過降低計算空間維度來減少存儲需求,在二維球坐標系下用FDTD計算閃電回擊通道產生的電場波形。MIGUEL等[12]從共享內存、分布式內存以及矢量化不同CPU體系結構等3個方面研究FDTD并行化計算,得出增加分布式節點,可提高計算性能等結論。MAI等用傳統FDTD計算粗網格區域的場值,用弱條件穩定時域有限差分(Weakly Conditionally Stable -FDTD)計算細網格,粗網格和細網格區域交界處采用了空間插值方法來減少仿真時間和內存需求[13],但過多的插值運算會降低計算精度。何欣波等[14]將Unconditionally Stable FDTD與傳統FDTD相結合克服了CFL(Courant-Friedrichs-Lewy)條件限制,由粗網格決定時間步長,在細網格區域及相鄰的粗網格區域計算特征值,通過減少特征值分解矩陣大小提高計算效率。DUAN等在電磁波的傳播方向上設置一個可移動窗口,可移動窗口內又設置了一個小窗口用來存儲大部分電磁場分量,當小窗口到達移動窗口邊界時把小窗口數據復制到新的移動窗口開始新區域的計算[15]。AHDAB等[16]為降低計算資源,把全部的障礙物設在激勵源附近區域,用FDTD計算出該區域內電磁波相互影響的全部結果,遠場電磁波則用拋物線方程計算,BAKIRTZIS等[17]則提出了一種FDTD和射線追蹤相結合的電磁波傳播方法來降低FDTD的算力需求。
FDTD方法可以分析巷道中防爆電氣開關、機電設備、金屬管道、動力和通信線纜等對礦井電磁波傳輸影響。但采用FDTD計算數百米三維立體巷道電磁波傳輸衰減,對計算設備的內存、算力都提出了非常高要求,制約了FDTD在礦井電磁波衰減仿真的應用。針對煤礦井下巷道斷面小,軸向長的特點,筆者提出了單位窗口循環使用的時域有限差分方法(Unit Window Recycling,UWR-FDTD)。該方法沿巷道軸向把整個巷道均勻劃分成多個大小相同的虛擬窗口,僅為一個與虛擬窗口大小相同的單位窗口分配內存空間,所有虛擬窗口根據自身編號循環使用單位窗口內存進行FDTD迭代,每一個虛擬窗口迭代完成后,存儲單位窗口電磁場,用于合成整個巷道內的電磁場。UWR-FDTD方法不改變FDTD的迭代方程,可根據計算機硬件設定單位窗口大小,劃分窗口個數,FDTD迭代過程僅在一個單位窗口的內存空間內進行,大幅減少了傳統FDTD方法對內存的需求,提高了計算效率。
麥克斯韋方程組揭示了電和磁之間的轉換規律,預測了電磁波的存在。FDTD直接求解麥克斯韋方程組研究電和磁之間的關系,成為計算電磁學中的一個重要數值計算方法[18]。FDTD把計算空間劃分成多個Yee網格,電場可以設置在Yee網格棱的中心,磁場在Yee網格面的中心,這樣每個電場分量周圍都環繞4個磁場分量,每個磁場分量也環繞著4個電場分量,與自然界中電場與磁場的物理結構保持一致。由于電場與磁場在時間和空間上都交叉分布,迭代過程也很容易用差分格式實現。用麥克斯韋方程組的兩個旋度方程,在時間和空間上差分可以推導出FDTD的迭代方程[19]。
(1)

式(1)中的E和H在直角坐標系下沿x,y,z三個方向可以拆分成關于Ex,Ey,Ez和Hx,Hy,Hz的6個標量方程。設定空間離散步長和時間離散步長,把上述6個標量方程中關于x,y,z,t的偏導數用中心差分近似后,即可得到電磁場在x,y,z三個方向分量上FDTD顯示表達式。FDTD的計算空間包括待仿真的全部空間以及空間內的激勵源和散射體,空間離散步長一般小于待研究頻域內最短波長的1/10[16]。受計算機內存容量限制,在模擬無限大區域內的電磁特性時,還需要在相應的截斷邊界處給出吸收邊界條件[21]。
對于斷面為拱形或馬蹄形的礦井巷道可以用面積等效為一個矩形巷道,故在直角坐標系下設等效矩形巷道斷面寬和高的大小分別為lx,ly,巷道軸向長lz,3個方向Yee網格空間步長均為δ。UWR-FDTD仿真時,用激勵源模擬發射天線,可以把電磁波信號通過更新方程耦合到空間網格內。激勵源所在的巷道斷面為軸向起始參考位置,一般情況下激勵源位于巷道斷面中心,激勵源在巷道斷面上的不同位置可以研究發射天線位置對礦井電磁波傳播特性的影響。沿巷道軸向把巷道劃分成多個長度大小均為lwz并從1順序編號的虛擬窗口,虛擬窗口的斷面大小與巷道的斷面大小相同。單位窗口中lwz的大小可根據計算機硬件性能和巷道軸向長度lz的大小來設定,一旦確定lwz后對計算結果向上取整,得單位窗口的數量nw:
(2)
仿真時,整個計算空間只需為電場和磁場分量及中間變量分配一個單位窗口的存儲空間用于FDTD迭代過程、與電磁波傳播lwz距離所用時間Tu相關的緩存空間。實際計算的巷道空間在3個方向上的大小分別為lx,ly和lz×nw,整個FDTD迭代過程在一個單位窗口空間內完成。電磁波在一個單位窗口內沿z方向的傳播時間Tu為
(3)
式中,Δt為Yee網格的時間步長。
電磁波從激勵源經過Tu時間到達單位窗口z方向的最后一個Yee網格,在此之后需要存儲最后一個網格的電場和磁場分量作為下一個虛擬窗口迭代的激勵源。
針對三維礦井巷道特點提出的UWR-FDTD方法同樣適用于一維和二維空間。使用UWR-FDTD方法對一維空間仿真時,Yee網格可比作一個點的擴展,需要按照時間步存儲最后一個點的各電場和磁場分量作為下一窗口迭代的子激勵源。若使用UWR-FDTD方法對二維空間進行仿真時,Yee網格可比作一條線的擴展,需要按照時間步存儲最后一條線的各電場和磁場分量作為下一窗口迭代的子激勵源。同理使用UWR-FDTD對三維空間進行仿真,Yee網格可作為一個面的擴展,需要按照時間步存儲最后一個面的各電場和磁場分量作為下一窗口迭代的子激勵源。電磁波在空氣介質中以近似光速從激勵源到達第一個窗口最后一個Yee網格的時間為Tu,并且激勵源以固定頻率周期變化,計算空間內各位置的電磁場也將呈周期性變化,在此可緩存單位窗口內最后一個網格Tu/Δt個時間步的電磁場作為下一個虛擬窗口的激勵源。每一個單位窗口運行完最后一個時間步長后都把該窗口的電場各分量轉存到合成電場中。UWR-FDTD的基本原理如圖1所示。

圖1 UWR-FDTD原理Fig.1 Principle of UWR-FDTD
為避免FDTD多個時間步迭代后的仿真結果發散,空間離散步長和時間離散步長上需要滿足CFL穩定性條件[19-20]:
(4)
若網格3個方向上的空間步長均為δ,則式(4)可簡化為
(5)
S=cΔt/δ
(6)
在直角坐標系中用x,y,z三個方向電場矢量和(下文用x方向表示計算空間的寬度方向,y方向表示高度方向,z方向表示長度方向)表示Yee網格的合成電場,該電場選在Yee網格的左下角。在兩個虛擬窗口或傳播介質與PML吸收邊界的交界面上,若使用左下角合成電場表示邊界處的電場,用最后一個Yee網格的電場表示傳播介質和PML交界處的電場,則存在一個Yee網格空間步長的計算誤差,如圖2中用Em(1,1,nm)表示第nm個網格的電場情況。因此,需要對邊界處的電場進一步處理,減小該誤差。邊界電場處理方法如圖2所示。

圖2 邊界電場處理方法Fig.2 Treatment method of boundary electric fields
圖2中傳播介質空間和PML空間沿x,y方向都包含nxny個Yee網格,傳播介質空間沿z方向分別被劃分為nm個網格,PML包含na個網格。傳播介質空間和PML空間第1個Yee網格的電場分別為Em(1,1,1)和Ea(1,1,1),傳播介質中最后一個Yee網格電場為Em(1,1,nm)。可以近似求解Em(1,1,nm+1)的值當作傳播介質和PML交界面處的電場,以及進入PML層的初始電場Ea(1,1,1),實現減少1個空間網格步長的誤差。FDTD迭代過程可得在兩種介質交界面上x和y方向上的電場Emx(1,1,nm+1)和Emy(1,1,nm+1),由于Emz(1,1,nm+1)在PML空間無法用傳播介質空間的電磁場分量直接計算,并且電磁波進入到PML空間后電磁場強度的衰減速率按照指數增加,在此采用傳播介質中Emz(1,1,nm-1)和Emz(1,1,nm)的線性插值表示Emz(1,1,nm+1)。可得在交界面上沿坐標軸3個方向上的電場為
(7)
本文實驗仿真軟件為Matlab R2021a,操作系統為Windows 10 專業版,CPU型號為Intel(R) Core(TM) i7-7700K CPU @ 4.20 GHz,系統內存容量32 GB。在此設計不同的實驗場景從計算精度與計算效率2個方面對比UWR-FDTD和傳統FDTD的數值計算性能,并研究UWR-FDTD中窗口大小與仿真時間之間的關系。
煤礦巷道壁材料主要為混凝土,其電磁參數變化不大,可以等效為各向同性的一個狹長的計算空間。設巷道斷面大小為nxny個Yee網格,軸向長為nz個Yee網格,所需電磁場各分量及中間變量個數為ne,則傳統FDTD所需內存空間κFDTD為
κFDTD=nxnynzne
(8)

(9)
單位窗口內循環計算出的電場分量最后要合成整個巷道的電場,因此還需要分配一個額外內存空間存儲整個巷道內的電場κa,大小為
κa=nxnynz
(10)
由式(9)和(10)可得UWR-FDTD所需的全部內存空間κUWR為
(11)
式(8)給出了傳統FDTD所需的內存空間,用式(11)除以式(8)可得UWR-FDTD所需的內存空間與傳統FDTD所需內存空間的比值,經過化簡后為
(12)
由式(12)可知ζ的大小與劃分窗口的個數以及電磁場及中間變量的個數有關。
下面在一維二維三維空間內分別分析。根據FDTD更新方程,一維計算空間需1個電場和1個磁場變量,用指數差分BerengerPML吸收邊界條件至少需4個系數變量,故ne≥6。二維計算空間需2個電場、2個磁場變量、至少4個參與FDTD迭代的系數變量,故二維計算空間中ne≥8。對于三維計算空間,電磁場將分裂成12個電磁場變量,以及至少4個系數變量,故三維計算空間中ne≥16。不同維度空間所需的變量個數不同,劃分不同窗口數量時UWR-FDTD方法和傳統FDTD方法所需內存的比值見表1。

表1 不同維度空間下ζ的大小
從表1可知,與傳統FDTD方法相比,用UWR-FDTD方法劃分的窗口數目越多,仿真時所需的內存空間越小;在三維計算空間內使用UWR-FDTD方法求解電磁問題時,節省的內存空間遠大于一維和二維計算空間。
FDTD仿真過程需要用與時間相關的函數作為激勵源引入到計算區域,常用周期時諧函數和脈沖函數作為時諧源和脈沖源。時諧源帶寬窄,一次可以仿真一個頻點的電磁波,但可以直觀的看出電磁波的傳播特性。脈沖源具有寬頻帶特性,使用脈沖源進行一次FDTD仿真,可以得到較寬頻帶上的電磁響應特性[22]。在此采用正弦波作為時諧源,微分高斯脈沖作為脈沖源對UWR-FDTD方法的性能進行分析,2種激勵源時間變化的函數方程為
(13)

對于式(13)中的微分高斯脈沖,當t0=τ且fmaxτ=2(fmax為微分高斯脈沖中待分析的最高電磁波頻率),可以用穩定因子S和單位波長包含的Yee網格數量dλ把Δt隱式表示,微分高斯脈沖可以離散為
(14)
在一維空間內待求計算空間長度為5 m,窗口長度為2 m,根據UWR-FDTD方法整個計算空間將被劃分為3個窗口,實際計算空間長度為6 m。Yee網格空間步長為載波波長的1/20[23],CFL穩定因子S為0.5,采用10層指數差分Berenger PML吸收邊界條件。以激勵源位置為起始點,用900 MHz的正弦波激勵源和微分高斯脈沖激勵,在無耗介質中迭代720時間步后電場的時域波形如圖3所示。由于迭代720時間步電場剛傳播到整個計算空間的最后一個網格,正弦波激勵源最右側電場還沒有完全穩定,使用UWR-FDTD方法劃分3個窗口合并后的電場與傳統的FDTD方法計算出的電場一致性較好,在2,4 m等2個窗口的交界處合成電場較為平滑。同一計算空間微分高斯脈沖激勵源在無耗介質中用UWR-FDTD和FDTD方法計算出的時域電場波形也幾乎完全一致。

圖3 一維無耗介質中的時域電場Fig.3 Time domain electric field in one-dimensionallossless medium
同無耗介質計算空間大小一樣,設有耗介質電導率σ為0.36 S/m,900 MHz正弦波和微分高斯脈沖激勵源傳播6 m,劃分3個窗口的UWR-FDTD與FDTD方法計算出的時域電場波形也一致,如圖4所示。

圖4 一維有耗介質中的時域電場Fig.4 Time domain electric field in one-dimensional loss medium
在長10 m、寬4 m的二維無耗介質計算空間中,設置UWR-FDTD單位窗口長3 m,寬4 m,可知需要把計算空間劃分4個虛擬窗口,實際計算空間長12 m、寬4 m。設正弦激勵源頻率為900 MHz,Yee網格空間步長16.67 mm,CFL穩定因子為0.5,時間步長27.8 ps,窗口四周采用10層指數差分BerengerPML吸收邊界條件。激勵源位于計算空間寬度方向上的中心位置,從激勵源開始在整個二維無耗計算空間內沿計算空間長度方向的時域電場分布如圖5所示。圖5(a),(b)為2種方法求解出在整個二維空間內的瞬時電場強度分布情況,圖5(c)為接收點位于寬度方向上的中心位置,沿巷道長度方向上各空間采樣點處的瞬時電場強度。

圖5 二維無耗介質中的時域電場Fig.5 Time domain electric field in two-dimensionallossless medium
由圖5可以看出,UWR-FDTD方法劃分4個窗口單獨運行后再合成的電場分布,與在一個窗口內使用FDTD運行的電場分布,在幅值和周期兩個方面都能較好吻合,表明UWR-FDTD能夠在二維計算空間內取得較好計算結果。
用文獻[2]中高3.4 m、寬4.8 m的一個等效矩形巷道來仿真三維無耗計算空間,窗口四周采用10層指數差分BerengerPML吸收邊界條件。當正弦激勵源頻率為900 MHz時,波長為333.11 mm,CFL穩定因子為0.40時,時間步長為22.22 ps,Yee網格空間步長為電磁波波長的1/20,對應單位窗口寬高分別288和204個Yee網格空間步。把偶極子激勵源放在巷道斷面中心,設單位窗口長為300個Yee網格,巷道長度約為5 m,用傳統FDTD與UWR-FDTD兩種方法,對3個虛擬窗口組成的巷道和5個虛擬窗口組成的巷道進行仿真,接收位置同樣位于巷道斷面中心沿巷道軸向采樣,最終得到各采樣點處的合成電場時域波形如圖6所示。
由圖6可知,在三維巷道內,用UWR-FDTD與傳統FDTD方法計算出的電場大小和周期較為一致,劃分多個虛擬窗口分別計算,最后合成的波形在2個窗口過渡時較為平滑,無明顯窗口分割痕跡。經過在不同維度空間對UWR-FDTD方法計算精度的分析可知,該方法可用于煤礦井下、金屬和非金屬等其他礦井、鐵路隧道、公路隧道、地鐵、地下人防工程、城市地下管廊、地下商場等空間內的電磁波分析。

圖6 三維無耗介質中時域電場Fig.6 Time domain electric field in three-dimensionallossless medium
UWR-FDTD方法把計算空間分割成多個虛擬窗口,電磁場分量及中間變量矩陣大小和單位窗口大小相同。傳統FDTD在1個時間步內迭代整個計算空間的電磁場變量矩陣,兩者相比,UWR-FDTD在每個時間步的計算量大幅減小,理論上可以縮短仿真時間,提高計算效率,下面在三維巷道內用仿真實驗給予驗證。
《煤礦安全規程》第九十條和《煤礦巷道斷面和交岔點設計規范》4.1.1節規定了軌道機車運輸的巷道凈高自軌道道面起不得低于2 m,采(盤)區內上山、下山和平巷的凈高度不得低于2m,薄煤層內不得低于1.8 m,巷道凈寬不宜小于2.0 m。用900 MHz正弦波作為激勵源,穩定因子S取0.4,空間步長為電磁波波長的1/20,約為16.67 mm。設巷道斷面寬高分別為140,110個Yee網格,約為2.33 m和1.84 m,單位窗口長為140個Yee網格約為2.33 m。測試從1到22個單位窗口組成的三維巷道,長度從2.33 m到51.33 m,在上述不同長度的巷道內用UWR-FDTD和FDTD兩種方法仿真同一電磁波傳播問題,仿真過程所需時間如圖7所示。

圖7 不同三維空間的仿真時間Fig.7 Simulation time in different three-dimensional spaces
當單位窗口長度和待仿真巷道長度相等時,相當于沒有劃分窗口,2種方法的仿真時間都為74 s,時間差為0,說明UWR-FDTD在不劃分窗口時與傳統FDTD的計算方法本質相同。當計算空間大小為8個虛擬窗口時,UWR-FDTD方法仿真用時2 739 s,FDTD方法仿真用時4 149 s,兩者用時相差1 410 s。當計算空間長度為22個虛擬窗口時,UWR-FDTD方法仿真用時19 499 s,FDTD方法仿真用時34 335 s,UWR-FDTD方法所需仿真時間比FDTD減少了14 836 s,降低了43.2%。故,UWR-FDTD方法非常適合狹長空間的電磁問題數值仿真,計算空間越大,同FDTD相比UWR-FDTD方法的計算效率越高。
對于同一段巷道,單位窗口的大小決定了需要劃分虛擬窗口的數量,虛擬窗口數量又決定著每個虛擬窗口在單位窗口內的迭代次數,迭代次數進一步影響整個巷道的仿真時間。因此需要研究巷道斷面的大小、巷道長度、電磁波頻率等多個條件下,虛擬窗口數量與仿真時間之間的關系,給出最優單位窗口大小的計算方法。
2.4.1 巷道斷面
用900 MHz正弦激勵源模擬發射天線,激勵源位于巷道斷面中心,Yee網格空間步長為波長的1/20,巷道長(z向)1 200個Yee網格,沿長度方向分別劃分為2,3,4,5,6,8,10,12,15,20,24個窗口,用xy表示巷道寬和高的Yee網格數量。當xy分別為140×110,90×70和50×40時,仿真上述11種虛擬窗口劃分方法所需的仿真時間見表2。

表2 不同巷道斷面大小的仿真時間
取表2中每一種xy仿真時間最小的一個虛擬窗口為時間基準(加粗數據),其他虛擬窗口的仿真時間都減去基準虛擬窗口的仿真時間,可得虛擬窗口數量、xy大小與仿真時間之間的關系,如圖8所示。
由圖8可以看出,對于同一長度的巷道,不管xy大小如何,劃分虛擬窗口數量過多或過少都會增加仿真時間,并且xy越大,劃分虛擬窗口的數量對仿真時間的影響越大。長度為1 200個Yee網格的巷道劃分8~10虛擬窗口時所用的仿真時間相對較少,此時單位窗口的長度為總巷道長度的1/10~1/8。
2.4.2 巷道長度
激勵源、穩定因子和吸收邊界條件與2.4.1節中設定的數值相同,巷道斷面固定為140×110個Yee網格,巷道長度分別為1 200,840和720個Yee網格。每種長度巷道依舊劃分為2,3,4,5,6,8,10,12,15,20,24個虛擬窗口,使用UWR-FDTD方法仿真上述不同長度巷道所需的仿真時間見表3。

表3 不同巷道長度的仿真時間
由表3可知,巷道長度方向Yee網格大小為1 200和840時,把上述巷道劃分為8個和5個虛擬窗口所需仿真時間最短;長度為720個Yee網格,劃分6個虛擬窗口時所需仿真時間最短。根據表3數據,以仿真時間最短的虛擬窗口為基準(加粗數據),其他虛擬窗口的仿真時間都減去基準窗口的仿真時間,可得不同巷道長度與仿真時間的關系,如圖9所示。

圖9 巷道長度與仿真時間的關系Fig.9 Relationship betweentunnel length and simulation time
由圖9可知不管待計算的巷道有多長,劃分虛擬窗口數量過多或過少都會增加仿真時間。對于同一xy斷面,單位窗口長度為對應巷道長度的1/10~1/6時,所需仿真用時相對較小。
2.4.3 電磁波頻率
用1 200,900,700 MHz的正弦波作為激勵源在長×寬×高為840×90×70個Yee網格的巷道內仿真,網格步長都為電磁波波長的1/10,計算不同電磁波頻率在劃分不同窗口數量時所需的仿真時間,統計劃分窗口數量與電磁波頻率之間的關系。3種電磁波頻率劃分不同窗口數量時所用的仿真時間見表4。

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

圖10 窗口數量和頻率與仿真時間的關系Fig.10 Relationship between electromagnetic wavefrequencies and simulation time
綜上所述,從巷道空間分析最優單位窗口的劃分方法,由圖8可知,不同巷道斷面單位窗口長度為巷道總長度的1/10~1/6內時,仿真時間最短;由圖9可知對于不同巷道長度,單位窗口長度也是在上述取值范圍內仿真時間最短。可得出在同一計算空間內,單位窗口長度為計算空間總長度的1/10~1/6時仿真計算用時最短。對于不同電磁波頻率,從圖10中可得當單位窗口長度為計算空間總長度的1/10~1/6時,計算效率較高。因此,在同一巷道內,不管電磁波頻率多大,當單位窗口長度為計算空間總長度的1/10~1/6時,都能取得較好的計算效率。
(1)在一維二維和三維計算空間內,不管是無耗介質還是有耗介質,UWR-FDTD計算出的數值結果與傳統FDTD計算出的數值結果都一致,具有較高的計算精度。
(2)同傳統FDTD方法相比,用UWR-FDTD方法劃分的窗口數目越多,仿真時所需的內存空間越小;在三維巷道中UWR-FDTD方法節省的內存空間遠大于一維和二維空間。
(3)同一巷道內UWR-FDTD所需的仿真時間與劃分虛擬窗口的數量、巷道斷面面積、巷道長度、電磁波頻率有關。單位窗口長度為待計算巷道總長度的1/10~1/6時,UWR-FDTD方法所需的仿真時間最短。