基于懸掛變量的顯式無條件穩定時域有限差分亞網格算法*

何欣波 魏兵

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

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

1 引言

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

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

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

2 EUS-FDTD 方法

對于無耗情形,電場的二階偏微分方程可寫為[1]

其中{e}表示計算區域中所有電場構成的向量,{f}為激勵源項,

為空間離散算子1/μ·(?×)·1/ε·(?×) 形成的系統矩陣,該矩陣只和空間網格離散尺寸以及網格的介質參數相關.其中,Se表示電場旋度離散后所構成的稀疏矩陣,Sh表示磁場旋度離散后所構成的稀疏矩陣.D1/ε表示元素為1/ε的對角矩陣,D1/μ表示元素為1/μ的對角矩陣.

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

一般而言,系統矩陣S中存在著穩定模式以及不穩定模式.穩定模式對應的特征值滿足

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

從本質上講,傳統FDTD 方法中的CFL 條件保證了數值系統中的所有模式均是穩定的.通常,傳統FDTD 方法在對目標進行網格離散時,僅需1/10 波長的網格離散尺度即可滿足精度要求.然而,當計算區域中含有精細結構時,需要選取較小的網格離散尺寸以便精確模擬精細結構的幾何形狀.由于CFL 條件的限制,需要選取較小的時間步長保證FDTD 算法時間迭代的穩定性.此時,小時間步長所對應的最高求解頻率遠大于所關心的最高求解頻率,造成了計算資源的浪費.

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

1)網格離散,構造系統矩陣S;

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

3)根據求解的最高頻率確定時間步長,找出不滿足式(4)的特征值以及對應的特征模式Vh;

4)濾除系統矩陣S中的不穩定模式,得到僅含穩定模式的系統矩陣Sl,其中

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

5)步驟4)保證了數值系統中僅含有穩定模式,在這些穩定模式中還存在著零空間模式,這些零空間模式會影響計算結果的正確性,因此還需要濾除零空間模式[17]:

通過以上步驟,即可實現大時間步下FDTD算法的穩定性.

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

3 懸掛變量亞網格技術

3.1 實施過程

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

圖1 懸掛變量亞網格技術Fig.1.The subgridding technique of hanging variables.

對于邊界處的粗網格電場,有

類似地,對于邊界處的細網格電場,有

對粗細網格邊界處的場采用如下的插值方式[19,20]:

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

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

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

對(14)式化簡可得

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

當Ratio=2 時,(15)式可簡化為

此外,邊界處細網格的磁場可寫為

(16)式和(20)式表示了邊界處電場和磁場之間的耦合.在其余區域,采用和均勻網格類似的方式構造系統矩陣.

3.2 穩定性證明

根據文獻[15]可知,若EUS-FDTD 方法所構成的系統矩陣是對稱半正定的,則系統矩陣的特征值為非負實數,此時我們可以根據所選擇的時間步長,濾除系統矩陣中大特征值對應的特征模式,實現無條件穩定.需要強調的是,對于亞網格數值系統,時間步長可由粗網格尺寸根據CFL 條件確定,此時EUSFDTD 算法具有較高的效率.本節從系統矩陣的對稱性角度出發,證明懸掛變量亞網格技術的穩定性.

為了簡單起見,對圖1 中的電磁場重新編號并消除冗余編號,得到如圖2 所示的電磁場編號.

圖2 編號重排后的懸掛變量亞網格技術Fig.2.The subgridding technique of hanging variables after numbering rearrangement.

因此,粗-細網格界面附近的磁場可寫為

粗-細網格界面附近的電場為

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

因此,最終所構成的系統矩陣為非對稱的.接下來,將給出系統矩陣的對稱化方法.

(29)式可寫為

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

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

因此,新的數值系統可寫為

4 顯式無條件穩定迭代

根據第3 節的推導,此時,關于電場的二階偏微分方程可寫為

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

5 數值算例

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

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

表1 三種數值算法的計算時間比較Table 1.Comparison of calculation times of three numerical algorithms.

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

圖4 直線AB 上采樣點的磁場強度Fig.4.Magnetic field intensity of sampling points on line AB.

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

5.2 多個介質目標的電磁散射特性

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

表2 三種數值算法的計算參數比較Table 2.Comparison of calculation parameters of three numerical algorithms.

圖6 多個介質目標計算模型Fig.6.Calculation model for multi medium targets.

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

圖8 監測點時域波形(a) P1;(b) P2Fig.8.Time domain waveform of monitoring points:(a) P1;(b) P2.

圖9 大時間步情形Uniform-FDTD 與HVS-FDTD 算法的時域波形(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 三維含介質PEC 腔體

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

表3 三種數值算法的計算參數比較Table 3.Comparison of calculation parameters of three numerical algorithms.

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

圖11 監測點時域波形Fig.11.Time domain waveform of monitoring the point.

6 結 論

本文提出了基于懸掛變量的EUS-FDTD 亞網格算法用于高效高精度分析含精細結構的電磁特性.本文給出了懸掛變量亞網格算法在邊界處的場值交換過程,從系統矩陣對稱性出發證明了懸掛變量亞網格算法的穩定性,并通過濾除亞網格數值系統中的不穩定模式實現顯式無條件穩定迭代.該方法具有實施過程簡單、精度高、效率高的優勢.本文所給的三個數值算例表明了所提方法的有效性.