基于懸掛變量的顯式無條件穩(wěn)定時(shí)域有限差分亞網(wǎng)格算法*

何欣波 魏兵

(西安電子科技大學(xué)物理學(xué)院,西安 710071)

時(shí)域有限差分(finite-difference time-domain,FDTD)方法由于穩(wěn)定性條件的限制,在處理含精細(xì)結(jié)構(gòu)的電磁問題時(shí)效率不高.顯式無條件穩(wěn)定(explicit unconditionally stable,EUS) FDTD 方法通過濾除系統(tǒng)矩陣的不穩(wěn)定模式,能夠消除穩(wěn)定性條件的限制,提高精細(xì)結(jié)構(gòu)的仿真效率.然而,EUS-FDTD 方法需要求解數(shù)值系統(tǒng)矩陣的特征值,在采用亞網(wǎng)格方案對(duì)含精細(xì)結(jié)構(gòu)目標(biāo)離散時(shí)需要保證數(shù)值系統(tǒng)矩陣的對(duì)稱性.現(xiàn)有的EUS-FDTD 亞網(wǎng)格方法存在著構(gòu)造復(fù)雜、精度不足等問題.針對(duì)以上問題,本文將懸掛變量亞網(wǎng)格(hanging variables subgridding,HVS)算法應(yīng)用于EUS-FDTD 算法中,從系統(tǒng)矩陣的對(duì)稱性出發(fā)證明了懸掛變量亞網(wǎng)格算法的穩(wěn)定性,給出了高精度、穩(wěn)定、容易實(shí)施的HVS-EUS-FDTD 方案.通過線磁流在自由空間的輻射、多個(gè)介質(zhì)目標(biāo)以及三維含介質(zhì)腔體的數(shù)值算例證明了所提方法的穩(wěn)定性、高精度以及高效性.數(shù)值實(shí)驗(yàn)表明,HVS-EUS-FDTD 算法的計(jì)算效率相比于均勻細(xì)網(wǎng)格FDTD 算法可提升數(shù)百倍,相比于HVS-FDTD 算法最高可提升Ratio (粗細(xì)網(wǎng)格尺寸比)倍.

1 引言

時(shí)域有限差分(finite-difference time-domain,FDTD)方法具有簡(jiǎn)單、高效、能夠處理復(fù)雜介質(zhì)的特點(diǎn),在計(jì)算電磁學(xué)中廣為流行[1-3].然而,FDTD 方法采用顯式迭代方案,其迭代時(shí)間步長(zhǎng)無法突破穩(wěn)定性條件Courant-Friedrichs-Lewy(CFL)的限制.在用FDTD 方法處理含電小尺寸結(jié)構(gòu)的電磁問題時(shí),需要使用非常小的網(wǎng)格對(duì)其進(jìn)行網(wǎng)格剖分,小網(wǎng)格尺寸將會(huì)造成小的時(shí)間迭代步長(zhǎng),進(jìn)而導(dǎo)致迭代步數(shù)巨大,仿真效率低下.為了消除或減弱穩(wěn)定性條件的限制,交替方向隱式(alternating direction implicit,ADI)[4-6]、局部一維(locally one-dimensional,LOD)[7-9]、克蘭克-尼科爾森(Crank-Nicolson,CN)[10,11]、弱條件穩(wěn)定(weakly conditionally stable,WCS)[12,13]等改進(jìn)的FDTD 方法被相繼提出,在某種程度上提升了FDTD 方法仿真含精細(xì)結(jié)構(gòu)目標(biāo)的計(jì)算效率.然而,上述方法要么需要求解矩陣方程,要么數(shù)值色散誤差較大,在實(shí)際應(yīng)用中需要平衡計(jì)算效率和精度.近年來,顯式無條件穩(wěn)定(explicit unconditionally stable,EUS) FDTD 方法被提出[14-17],該方法通過濾除系統(tǒng)矩陣中大特征值對(duì)應(yīng)的特征模式實(shí)現(xiàn)無條件穩(wěn)定,保留了傳統(tǒng)FDTD 方法顯式迭代的優(yōu)點(diǎn),并且具有較高的計(jì)算精度.

EUS-FDTD 方法和傳統(tǒng)FDTD 方法類似,采用結(jié)構(gòu)網(wǎng)格(Yee 元胞)離散計(jì)算區(qū)域.當(dāng)計(jì)算區(qū)域中含有電小尺寸目標(biāo)時(shí),為了精確模擬該目標(biāo)的幾何結(jié)構(gòu),需要使用非常小的網(wǎng)格對(duì)其進(jìn)行離散.一種較為簡(jiǎn)單粗暴的方式是對(duì)整個(gè)計(jì)算區(qū)域采用均勻細(xì)網(wǎng)格離散,這會(huì)造成計(jì)算區(qū)域的網(wǎng)格規(guī)模急劇增多,計(jì)算效率較低.亞網(wǎng)格技術(shù)是一種離散含電小尺寸目標(biāo)的較好方式,然而,應(yīng)用亞網(wǎng)格技術(shù)時(shí)需要謹(jǐn)慎處理粗-細(xì)網(wǎng)格交界處的場(chǎng),否則會(huì)造成計(jì)算的不穩(wěn)定性.2017 年,Yan 和Jiao[15]通過構(gòu)造對(duì)稱半正定(symmetric positive semidefinite,SPD)算子,實(shí)現(xiàn)了FDTD 亞網(wǎng)格數(shù)值系統(tǒng)的對(duì)稱性,保證了計(jì)算的穩(wěn)定性,但是該方法構(gòu)造對(duì)稱半正定算子的過程較為復(fù)雜,并且在三維情形僅考慮了細(xì)網(wǎng)格區(qū)域的場(chǎng)耦合到粗網(wǎng)格區(qū)域中,并未考慮粗網(wǎng)格區(qū)域的場(chǎng)耦合到細(xì)網(wǎng)格區(qū)域中,這會(huì)影響計(jì)算精度.2018 年,Yan 和Jiao[16]通過三角形內(nèi)插方法構(gòu)造非對(duì)稱的FDTD 亞網(wǎng)格數(shù)值系統(tǒng),并采用后向差分離散保證計(jì)算的穩(wěn)定性,此外,還利用Neumann 級(jí)數(shù)實(shí)現(xiàn)后向差分FDTD 方法的顯式迭代.然而,相比于中心差分離散,后向差分精度不足,并且Neumann 級(jí)數(shù)展開會(huì)消耗大量的計(jì)算資源.2020 年,Zeng 和Jiao[18]提出了一種系統(tǒng)的方法在空間和時(shí)間上構(gòu)造FDTD 亞網(wǎng)格SPD 算子,保證了計(jì)算的穩(wěn)定性.與Yan 提出的SPD 算子類似,該方法構(gòu)造SPD 算子較為復(fù)雜,并且需要在粗細(xì)網(wǎng)格區(qū)域設(shè)置不同的時(shí)間步長(zhǎng)保證計(jì)算的穩(wěn)定性.2017 年,Bekmambetova 等[19]將耗散系統(tǒng)理論應(yīng)用于FDTD 方法中,并根據(jù)該理論給出了一種穩(wěn)定的、易于實(shí)施的FDTD 亞網(wǎng)格方案,稱為懸掛變量法(hanging variables,HV).該方法通過在粗細(xì)網(wǎng)格邊界處引入懸掛變量,使得粗網(wǎng)格區(qū)域和細(xì)網(wǎng)格區(qū)域變?yōu)閮蓚€(gè)單獨(dú)的子系統(tǒng),這兩個(gè)子系統(tǒng)之間通過懸掛變量耦合在一起.

本文將懸掛變量亞網(wǎng)格(hanging variables subgridding,HVS)算法引入到EUS-FDTD 算法中,并從系統(tǒng)矩陣的對(duì)稱性方面證明懸掛變量亞網(wǎng)格算法的穩(wěn)定性,通過濾除系統(tǒng)矩陣中大特征對(duì)應(yīng)的不穩(wěn)定模式,使得整個(gè)計(jì)算區(qū)域采用均勻的大時(shí)間步迭代,提高含電小結(jié)構(gòu)目標(biāo)電磁仿真的計(jì)算效率.本文結(jié)構(gòu)如下: 第2 節(jié)介紹了EUS-FDTD 算法的基本原理和實(shí)現(xiàn)步驟;第3 節(jié)給出了懸掛變量亞網(wǎng)格算法的實(shí)施過程,并通過系統(tǒng)矩陣的對(duì)稱性證明了該方法的穩(wěn)定性;第4 節(jié)給出了懸掛變量亞網(wǎng)格算法的顯式無條件穩(wěn)定迭代的過程;第5 節(jié)通過線磁流在自由空間的輻射、多個(gè)介質(zhì)目標(biāo)以及三維含介質(zhì)腔體的電磁散射特性證明了所提方法的有效性;第6 節(jié)為本文的結(jié)論.

2 EUS-FDTD 方法

對(duì)于無耗情形,電場(chǎng)的二階偏微分方程可寫為[1]

其中{e}表示計(jì)算區(qū)域中所有電場(chǎng)構(gòu)成的向量,{f}為激勵(lì)源項(xiàng),

為空間離散算子1/μ·(?×)·1/ε·(?×) 形成的系統(tǒng)矩陣,該矩陣只和空間網(wǎng)格離散尺寸以及網(wǎng)格的介質(zhì)參數(shù)相關(guān).其中,Se表示電場(chǎng)旋度離散后所構(gòu)成的稀疏矩陣,Sh表示磁場(chǎng)旋度離散后所構(gòu)成的稀疏矩陣.D1/ε表示元素為1/ε的對(duì)角矩陣,D1/μ表示元素為1/μ的對(duì)角矩陣.

對(duì)(1)式采用中心差分離散,可得到FDTD的時(shí)域步進(jìn)形式為

一般而言,系統(tǒng)矩陣S中存在著穩(wěn)定模式以及不穩(wěn)定模式.穩(wěn)定模式對(duì)應(yīng)的特征值滿足

其中Δt為時(shí)間步長(zhǎng),λ為系統(tǒng)矩陣S的特征值.當(dāng)所取的時(shí)間步長(zhǎng)滿足(4)式時(shí),數(shù)值系統(tǒng)中僅存在穩(wěn)定模式,迭代過程是穩(wěn)定的.在網(wǎng)格的介質(zhì)參數(shù)不變的情況下,網(wǎng)格尺寸越小,特征值λ越大,為了保證迭代過程的穩(wěn)定性,需要取較小的時(shí)間步長(zhǎng)確保系統(tǒng)矩陣中僅含有穩(wěn)定模式,即時(shí)間步長(zhǎng)受限于空間網(wǎng)格尺寸.需要注意的是,(4)式和CFL 條件是等價(jià)的.在分析含電小尺寸目標(biāo)的電磁問題時(shí),這種限制會(huì)極大影響計(jì)算效率.

從本質(zhì)上講,傳統(tǒng)FDTD 方法中的CFL 條件保證了數(shù)值系統(tǒng)中的所有模式均是穩(wěn)定的.通常,傳統(tǒng)FDTD 方法在對(duì)目標(biāo)進(jìn)行網(wǎng)格離散時(shí),僅需1/10 波長(zhǎng)的網(wǎng)格離散尺度即可滿足精度要求.然而,當(dāng)計(jì)算區(qū)域中含有精細(xì)結(jié)構(gòu)時(shí),需要選取較小的網(wǎng)格離散尺寸以便精確模擬精細(xì)結(jié)構(gòu)的幾何形狀.由于CFL 條件的限制,需要選取較小的時(shí)間步長(zhǎng)保證FDTD 算法時(shí)間迭代的穩(wěn)定性.此時(shí),小時(shí)間步長(zhǎng)所對(duì)應(yīng)的最高求解頻率遠(yuǎn)大于所關(guān)心的最高求解頻率,造成了計(jì)算資源的浪費(fèi).

EUS-FDTD 方法在選取時(shí)間步長(zhǎng)上與傳統(tǒng)FDTD 方法有所區(qū)別.EUS-FDTD 方法由最高求解頻率并利用采樣定理確定時(shí)間步長(zhǎng),通過求解系統(tǒng)矩陣的特征值和特征模式并將大特征值(不滿足(4)式)對(duì)應(yīng)的特征模式濾除,使得整個(gè)數(shù)值系統(tǒng)中僅含有穩(wěn)定模式,保證了大時(shí)間步的穩(wěn)定迭代.具體步驟如下:

1)網(wǎng)格離散,構(gòu)造系統(tǒng)矩陣S;

2)求解系統(tǒng)矩陣的部分特征值和特征向量;

3)根據(jù)求解的最高頻率確定時(shí)間步長(zhǎng),找出不滿足式(4)的特征值以及對(duì)應(yīng)的特征模式Vh;

4)濾除系統(tǒng)矩陣S中的不穩(wěn)定模式,得到僅含穩(wěn)定模式的系統(tǒng)矩陣Sl,其中

VhT表示Vh的轉(zhuǎn)置,在時(shí)域步進(jìn)過程中,需要將(3)式中的S替換為Sl;

5)步驟4)保證了數(shù)值系統(tǒng)中僅含有穩(wěn)定模式,在這些穩(wěn)定模式中還存在著零空間模式,這些零空間模式會(huì)影響計(jì)算結(jié)果的正確性,因此還需要濾除零空間模式[17]:

通過以上步驟,即可實(shí)現(xiàn)大時(shí)間步下FDTD算法的穩(wěn)定性.

從理論上講,在滿足采樣精度的情形下,EUSFDTD 方法時(shí)間步長(zhǎng)的選取無限制.但是,由于該方法需要求解系統(tǒng)矩陣的部分特征值和特征向量,當(dāng)時(shí)間步長(zhǎng)選取得很大,那就意味著系統(tǒng)矩陣中的不穩(wěn)定模式很多,此時(shí)求解特征值和特征向量這一步驟就會(huì)消耗大量的計(jì)算資源.一般而言,該方法更適用于含精細(xì)結(jié)構(gòu)的電磁仿真中.精細(xì)結(jié)構(gòu)采用細(xì)網(wǎng)格離散,其余大部分區(qū)域采用粗網(wǎng)格離散,時(shí)間步長(zhǎng)僅由粗網(wǎng)格尺寸確定,此時(shí)不穩(wěn)定模式僅存在細(xì)網(wǎng)格及其周圍一層粗網(wǎng)格中,通過部分特征值求解技術(shù)可快速獲取系統(tǒng)矩陣的不穩(wěn)定模式,此時(shí)該方法的計(jì)算效率較高.

3 懸掛變量亞網(wǎng)格技術(shù)

3.1 實(shí)施過程

當(dāng)計(jì)算區(qū)域中含有電小尺寸目標(biāo)時(shí),采用亞網(wǎng)格技術(shù)對(duì)離散計(jì)算區(qū)域能夠極大降低網(wǎng)格量,進(jìn)而提高計(jì)算效率.此時(shí),在粗-細(xì)網(wǎng)格邊界處,由于兩種網(wǎng)格的采樣點(diǎn)不一致,計(jì)算可能用到不在采樣點(diǎn)的場(chǎng)值,因此需要采用合理的插值方案將這些點(diǎn)的場(chǎng)值用已知采樣點(diǎn)的場(chǎng)值表示.本文以二維TE 情形粗-細(xì)網(wǎng)格邊界為例,討論插值方法的思想.對(duì)于TE 波,電場(chǎng)位于面片四條棱邊的中心,磁場(chǎng)位于面片的中心.粗網(wǎng)格尺寸為Δc,細(xì)網(wǎng)格尺寸為Δf,Δc與Δf的關(guān)系為Δc=Ratio·Δf,其中Ratio為粗細(xì)網(wǎng)格尺寸比,圖1 為Ratio=2 情形.交界處附近粗網(wǎng)格的電場(chǎng)和磁場(chǎng)分別用E1和H1表示,交界處附近細(xì)網(wǎng)格的電場(chǎng)和磁場(chǎng)分別用ei(i=1-7)和hi(i=1,2) 表示.在邊界處,采用懸掛變量連接界面處粗網(wǎng)格和細(xì)網(wǎng)格中的場(chǎng).在粗網(wǎng)格邊界處,引入關(guān)于磁場(chǎng)的懸掛變量.類似地,在細(xì)網(wǎng)格邊界處,引入關(guān)于磁場(chǎng)的懸掛變量(i=1,2)[19].在推導(dǎo)交界面處場(chǎng)的最終迭代式的過程中,這些懸掛變量將被消去.

圖1 懸掛變量亞網(wǎng)格技術(shù)Fig.1.The subgridding technique of hanging variables.

對(duì)于邊界處的粗網(wǎng)格電場(chǎng),有

類似地,對(duì)于邊界處的細(xì)網(wǎng)格電場(chǎng),有

對(duì)粗細(xì)網(wǎng)格邊界處的場(chǎng)采用如下的插值方式[19,20]:

由(7)式和(10)式可得

對(duì)(8)式中的所有項(xiàng)求和,并將(9)式代入可得

將(13)式代入到(11)式中,可得

對(duì)(14)式化簡(jiǎn)可得

粗網(wǎng)格上的電場(chǎng)E1計(jì)算完成后,可通過(9)式直接得到細(xì)網(wǎng)格邊界處的電場(chǎng).

當(dāng)Ratio=2 時(shí),(15)式可簡(jiǎn)化為

此外,邊界處細(xì)網(wǎng)格的磁場(chǎng)可寫為

(16)式和(20)式表示了邊界處電場(chǎng)和磁場(chǎng)之間的耦合.在其余區(qū)域,采用和均勻網(wǎng)格類似的方式構(gòu)造系統(tǒng)矩陣.

3.2 穩(wěn)定性證明

根據(jù)文獻(xiàn)[15]可知,若EUS-FDTD 方法所構(gòu)成的系統(tǒng)矩陣是對(duì)稱半正定的,則系統(tǒng)矩陣的特征值為非負(fù)實(shí)數(shù),此時(shí)我們可以根據(jù)所選擇的時(shí)間步長(zhǎng),濾除系統(tǒng)矩陣中大特征值對(duì)應(yīng)的特征模式,實(shí)現(xiàn)無條件穩(wěn)定.需要強(qiáng)調(diào)的是,對(duì)于亞網(wǎng)格數(shù)值系統(tǒng),時(shí)間步長(zhǎng)可由粗網(wǎng)格尺寸根據(jù)CFL 條件確定,此時(shí)EUSFDTD 算法具有較高的效率.本節(jié)從系統(tǒng)矩陣的對(duì)稱性角度出發(fā),證明懸掛變量亞網(wǎng)格技術(shù)的穩(wěn)定性.

為了簡(jiǎn)單起見,對(duì)圖1 中的電磁場(chǎng)重新編號(hào)并消除冗余編號(hào),得到如圖2 所示的電磁場(chǎng)編號(hào).

圖2 編號(hào)重排后的懸掛變量亞網(wǎng)格技術(shù)Fig.2.The subgridding technique of hanging variables after numbering rearrangement.

因此,粗-細(xì)網(wǎng)格界面附近的磁場(chǎng)可寫為

粗-細(xì)網(wǎng)格界面附近的電場(chǎng)為

由(26)式和(29)式可知,

因此,最終所構(gòu)成的系統(tǒng)矩陣為非對(duì)稱的.接下來,將給出系統(tǒng)矩陣的對(duì)稱化方法.

(29)式可寫為

將其代入到(24)式和(27)式中,此時(shí),(26)式重寫為

由(35)式和(36)式可知,

因此,新的數(shù)值系統(tǒng)可寫為

4 顯式無條件穩(wěn)定迭代

根據(jù)第3 節(jié)的推導(dǎo),此時(shí),關(guān)于電場(chǎng)的二階偏微分方程可寫為

在每一步迭代完成后還需要濾除零空間模式,即

5 數(shù)值算例

5.1 線磁流在自由空間的輻射

二維計(jì)算區(qū)域的大小為6 m×6 m,區(qū)域外采用吸收邊界截?cái)?在計(jì)算區(qū)域中心處設(shè)置了一個(gè)頻率為f=0.3 GHz的線磁流源.該計(jì)算模型采用了基于懸掛變量的亞網(wǎng)格方法進(jìn)行離散,其中粗網(wǎng)格尺寸為Δc=0.05 m,線磁流源周圍0.15 m 的方形區(qū)域內(nèi)采用細(xì)網(wǎng)格離散,粗細(xì)網(wǎng)格比為Ratio=5,即細(xì)網(wǎng)格尺寸為寸為Δf=0.01 m,計(jì)算模型示意如圖3 所示.在本例中,采用了解析解(analytical)、均勻細(xì)網(wǎng)格的傳統(tǒng)FDTD 方法(uniform FDTD),HVS-FDTD 及HVS-EUS-FDTD 計(jì)算了直線AB上采樣點(diǎn)的磁場(chǎng)強(qiáng)度,以上四種方法的計(jì)算結(jié)果如圖4 所示.由圖4 可知,三種數(shù)值方法與解析解的結(jié)果吻合得很好.此外,表1 給出了三種數(shù)值方法的計(jì)算時(shí)間、時(shí)間步長(zhǎng)和內(nèi)存占用情況.uniform FDTD 和HVS-FDTD 需要采用較小的時(shí)間步長(zhǎng)進(jìn)行迭代以保證計(jì)算的穩(wěn)定性,而HVS-EUSFDTD 可采用較大的時(shí)間步長(zhǎng).uniform FDTD 方法需要2295.5 s 完成仿真,采用懸掛變量亞網(wǎng)格算法進(jìn)行離散后,HVS-FDTD 僅需要17.45 s 即可完成仿真.更進(jìn)一步,將懸掛變量亞網(wǎng)格算法引入到EUS-FDTD 算法,HVS-EUS-FDTD 方法僅需要5.35 s 即可完成整個(gè)求解過程,其中特征值求解時(shí)間為0.05 s,迭代時(shí)間5.30 s.這是因?yàn)橄啾扔趗niform FDTD 方法,HVS-FDTD 采用亞網(wǎng)格離散,極大的降低了網(wǎng)格規(guī)模,提高了計(jì)算效率.而HVS-EUS-FDTD 在HVS-FDTD 方法的基礎(chǔ)上,濾除系統(tǒng)矩陣中的不穩(wěn)定模式,使得整個(gè)迭代過程可以采用均勻的大時(shí)間步進(jìn)行仿真,進(jìn)一步提高了計(jì)算效率.從uniform FDTD 和HVS-FDTD 方法的內(nèi)存占用情況也可以看出,亞網(wǎng)格方案極大降低了內(nèi)存,相比于HVS-FDTD 方法,HVS-EUSFDTD 方法需要額外存儲(chǔ)不穩(wěn)定模式,因此HVSEUS-FDTD 方法的內(nèi)存占用量比HVS-FDTD 方法的稍大.圖5 還給出了HVS-EUS-FDTD 方法中大特征值的分布,在本例中特征值征值λ ＞4 的特征值有210 個(gè),通過將其對(duì)應(yīng)的特征模式濾除即可實(shí)現(xiàn)整個(gè)數(shù)值系統(tǒng)的無條件穩(wěn)定性.

表1 三種數(shù)值算法的計(jì)算時(shí)間比較Table 1.Comparison of calculation times of three numerical algorithms.

圖3 計(jì)算模型示意圖Fig.3.The schematic diagram of calculation model.

圖4 直線AB 上采樣點(diǎn)的磁場(chǎng)強(qiáng)度Fig.4.Magnetic field intensity of sampling points on line AB.

圖5 特征值分布Fig.5.The distribution of eigenvalues.

5.2 多個(gè)介質(zhì)目標(biāo)的電磁散射特性

本例計(jì)算了TE 情形多個(gè)無限長(zhǎng)介質(zhì)圓柱的散射問題.計(jì)算模型如圖6 所示.在計(jì)算區(qū)域中有9 個(gè)相對(duì)介電常數(shù)εr=25 的無限長(zhǎng)介質(zhì)圓柱,圓柱的直徑d=0.05 m,圓柱之間的距離均為b=0.1 m .入射源為線磁流源,源的形式為微分高斯脈沖,脈沖寬度τ=1 ns .在計(jì)算區(qū)域內(nèi)設(shè)置有兩個(gè)監(jiān)測(cè)點(diǎn),分別為P1和P2.加源點(diǎn)以及監(jiān)測(cè)點(diǎn)位置如圖6 所示.本例中,采用了uniform FDTD,HVS-FDTD以及HVS-EUS-FDTD 三種方法對(duì)該模型進(jìn)行仿真.由于計(jì)算區(qū)域中存在介質(zhì),為了提高建模精度,需要采用細(xì)網(wǎng)格對(duì)介質(zhì)區(qū)域進(jìn)行網(wǎng)格離散.對(duì)于uniform FDTD 方法,采用均勻細(xì)網(wǎng)格對(duì)該模型進(jìn)行離散,細(xì)網(wǎng)格尺寸為Δf=0.002 m ;對(duì)于HVSFDTD 以及HVS-EUS-FDTD 方法,采用懸掛變量亞網(wǎng)格方法對(duì)該模型進(jìn)行離散,對(duì)于單個(gè)介質(zhì)目標(biāo)的網(wǎng)格離散示意圖如圖7 所示.對(duì)于懸掛變量亞網(wǎng)格,粗網(wǎng)格尺寸Δc=0.01 m,細(xì)網(wǎng)格尺寸與uniform FDTD 方法的相同,為Δf=0.002 m .以上三種方法所獲得的監(jiān)測(cè)點(diǎn)P1和P2的時(shí)域波形如圖8 所示,可以看到,HVS-FDTD 以及HVSEUS-FDTD 方法的結(jié)果與uniform FDTD 方法的結(jié)果吻合得很好,證明了所提方法的精度.為了獲得穩(wěn)定的仿真結(jié)果,uniform FDTD 方法的時(shí)間步長(zhǎng)需為3.33×10-12s,完成整個(gè)仿真需花費(fèi)91.36 s;HVS-FDTD 方法的計(jì)算區(qū)域中有細(xì)網(wǎng)格的存在,其時(shí)間步長(zhǎng)和uniform FDTD 的相同,但亞網(wǎng)格技術(shù)降低了網(wǎng)格規(guī)模,使得HVS-FDTD 方法僅需2.19 s 即可完成整個(gè)仿真;HVS-EUS-FDTD 方法通過濾除HVS-FDTD 方法的不穩(wěn)定模式,保證了整個(gè)計(jì)算區(qū)域采用了統(tǒng)一的大時(shí)間步進(jìn)行迭代,減少了迭代步數(shù),使得僅需1.07 s 即可完成計(jì)算,其中特征值求解花費(fèi)0.29 s,時(shí)域迭代花費(fèi)0.78 s.圖9 給出了時(shí)間步長(zhǎng)擴(kuò)大2 倍情形uniform FDTD和HVS-FDTD 方法所計(jì)算的監(jiān)測(cè)點(diǎn)時(shí)域波形,可以看到,此時(shí)uniform FDTD 和HVS-FDTD 方法不穩(wěn)定.在內(nèi)存消耗方面,uniform FDTD 方法消耗的內(nèi)存最多,HVS-FDTD 方法消耗的最少,HVSEUS-FDTD 方法消耗的內(nèi)存比HVS-FDTD方法稍大.三種方法的詳細(xì)計(jì)算參數(shù)如表2 所列.從以上可知,HVS-EUS-FDTD 相比于uniform FDTD以及HVS-FDTD 方法有更高的計(jì)算效率.

表2 三種數(shù)值算法的計(jì)算參數(shù)比較Table 2.Comparison of calculation parameters of three numerical algorithms.

圖6 多個(gè)介質(zhì)目標(biāo)計(jì)算模型Fig.6.Calculation model for multi medium targets.

圖7 懸掛變量亞網(wǎng)格離散模型Fig.7.Subgridding discrete model with hanging variables.

圖8 監(jiān)測(cè)點(diǎn)時(shí)域波形(a) P1;(b) P2Fig.8.Time domain waveform of monitoring points:(a) P1;(b) P2.

圖9 大時(shí)間步情形Uniform-FDTD 與HVS-FDTD 算法的時(shí)域波形(a)Uniform-FDTD;(b) HVS-FDTDFig.9.Time domain waveforms of Uniform-FDTD and HVS-FDTD algorithms with large time steps: (a)Uniform-FDTD;(b) HVSFDTD.

5.3 三維含介質(zhì)PEC 腔體

三維PEC 腔體尺寸為2 m×2 m×2 m,腔體中心區(qū)域含一塊大小為0.05 m×0.05 m×0.05 m、介質(zhì)參數(shù)為εr=20 的介質(zhì)塊,計(jì)算模型如圖10 所示.在PEC 腔體內(nèi)部(0.2 m,0.2 m,0.225 m)處設(shè)置z方向的電偶極子源,其形式為微分高斯脈沖,脈沖寬度τ=4 ns .在(0.7 m,0.7 m,0.725 m)處設(shè)置一個(gè)監(jiān)測(cè)點(diǎn).本例也采用了uniform FDTD,HVS-FDTD 以及HVS-EUS-FDTD 三種方法對(duì)該模型進(jìn)行仿真.對(duì)于uniform FDTD 方法,采用均勻細(xì)網(wǎng)格對(duì)該模型進(jìn)行離散,細(xì)網(wǎng)格尺寸為Δf=0.01 m;對(duì)于HVS-FDTD 以及HVS-EUSFDTD 方法,采用懸掛變量亞網(wǎng)格方法對(duì)該模型進(jìn)行離散,介質(zhì)塊為細(xì)網(wǎng)格,其余區(qū)域?yàn)榇志W(wǎng)格,其中粗網(wǎng)格尺寸為Δc=0.05 m,細(xì)網(wǎng)格尺寸與uniform FDTD 方法的相同,即Δf=0.01 m .以上三種方法所獲得的監(jiān)測(cè)點(diǎn)的時(shí)域波形如圖11 所示.可以看到,三種方法的結(jié)果符合得很好,證明了所提方法的精度.為了獲得穩(wěn)定的仿真結(jié)果,uniform FDTD 方法的時(shí)間步長(zhǎng)需為1.67×10-11s,完成整個(gè)仿真需花費(fèi)28816.8 s;HVS-FDTD 方法的計(jì)算區(qū)域中有細(xì)網(wǎng)格的存在,其時(shí)間步長(zhǎng)和uniform FDTD 的相同,但亞網(wǎng)格技術(shù)降低了網(wǎng)格規(guī)模,使得HVS-FDTD 方法僅需314.2 s 即可完成整個(gè)仿真;HVS-EUS-FDTD 方法通過濾除HVS-FDTD方法的不穩(wěn)定模式,保證了整個(gè)計(jì)算區(qū)域采用了統(tǒng)一的大時(shí)間步進(jìn)行迭代,減少了迭代步數(shù),使得僅需64.71 s 即可完成計(jì)算.uniform FDTD 方法消耗的內(nèi)存最多,HVS-FDTD 方法消耗的最少,HVSEUS-FDTD 方法消耗的內(nèi)存比HVS-FDTD 方法稍大.三種方法的詳細(xì)計(jì)算參數(shù)如表3 所列.從以上可知,在三維情形,HVS-EUS-FDTD 相比于uniform FDTD 以及HVS-FDTD的計(jì)算效率提升更明顯.

表3 三種數(shù)值算法的計(jì)算參數(shù)比較Table 3.Comparison of calculation parameters of three numerical algorithms.

圖10 三維含介質(zhì)PEC 腔體計(jì)算模型Fig.10.The computational model for 3-D PEC cavity containing a dielectric.

圖11 監(jiān)測(cè)點(diǎn)時(shí)域波形Fig.11.Time domain waveform of monitoring the point.

6 結(jié) 論

本文提出了基于懸掛變量的EUS-FDTD 亞網(wǎng)格算法用于高效高精度分析含精細(xì)結(jié)構(gòu)的電磁特性.本文給出了懸掛變量亞網(wǎng)格算法在邊界處的場(chǎng)值交換過程,從系統(tǒng)矩陣對(duì)稱性出發(fā)證明了懸掛變量亞網(wǎng)格算法的穩(wěn)定性,并通過濾除亞網(wǎng)格數(shù)值系統(tǒng)中的不穩(wěn)定模式實(shí)現(xiàn)顯式無條件穩(wěn)定迭代.該方法具有實(shí)施過程簡(jiǎn)單、精度高、效率高的優(yōu)勢(shì).本文所給的三個(gè)數(shù)值算例表明了所提方法的有效性.