李磊 衛(wèi)敏

（安徽大學電氣工程與自動化學院，合肥 230601）

引 言

時域多分辨率(multiresolution time-domain, MRTD)算法是將矩量法(moment of method, MoM) 中的小波伽遼金法和多分辨率分析(multiresolution analysis,MRA)相結(jié)合[1-2]，電磁場量在空間方向上用尺度函數(shù)和小波函數(shù)展開處理[3-5]，在時間方向上用矩形脈沖函數(shù)進行展開. 與時域有限差分(finitedifference time-domain, FDTD)算法相比，具有良好的色散特性，從而可節(jié)省計算資源，提高計算效率.但是穩(wěn)定度較低，對時間步的限制比FDTD算法更加的嚴苛，也就是“用時間換空間”. 文獻[6]提出的交替方向隱式(alternating-direction implicit,ADI)-MRTD算法擺脫了傳統(tǒng)MRTD對于穩(wěn)定性條件的束縛[6]，但是其數(shù)值色散特性較差. 龍格-庫塔(Runge-Kutta, RK)-MRTD算法[7]在保持MRTD算法較好色散特性的同時進一步提高了數(shù)值穩(wěn)定性，提高了計算效率，但該算法是能量耗散的. 本文把辛算法引入傳統(tǒng)的MRTD算法中，在空間方向利用小波尺度函數(shù)進行離散處理，以獲得良好的色散特性，在時間方向利用辛算子進行離散，從而改善MRTD算法對于穩(wěn)定性條件的限制.

波導結(jié)構(gòu)元器件在微波毫米波通信中具有廣泛的應(yīng)用，多年來波導結(jié)構(gòu)的特性分析一直是電磁工程的重要研究方向. 波導結(jié)構(gòu)主要包括波導內(nèi)部的不連續(xù)(類諧振腔、填充介質(zhì)、膜片、導體柱等)，矩形波導縫隙，脊波導及其縫隙等更加復(fù)雜的結(jié)構(gòu)[8]. 矩形波導是波導結(jié)構(gòu)元器件中的典型代表，因此，辛時域多分辨率(symplectic multiresolution time-domain, S-MRTD)算法在矩形波導仿真中的研究，對于分析和理解其他波導結(jié)構(gòu)特性時具有參考意義.

1 S-MRTD算法

1.1 Maxwell方程的Hamilton表示

式中：E是電場強度；H是磁場強度；ε是媒質(zhì)中的介電常數(shù)；μ是媒質(zhì)中磁導率.

將式(1)寫成矩陣形式：

式中： ?為3維旋度算子表示形式；算符U和V含有三維旋度算子 ?.

1.2 空間離散

在空間方向上對上述方程進行離散，采用MRTD算法將電場與磁場用具有緊支撐的Daubechies小波尺度函數(shù)進行展開并代入Maxwell方程組，得到

式中：LS為 基函數(shù)的支撐域；a(v)為Daubechies小波尺度函數(shù)的連接系數(shù).

1.3 時間離散

在時間方向上，利用辛算法進行離散處理，時間t取值為0～ ?t時，式(2)的解可表示為

由于指數(shù)算符 exp(?t(U+V))并不能顯式計算，式(7) 僅僅是形式解，因此，采用辛積分傳播子[9]得到 exp(?t(U+V))的m級p階近似展開為

式中：m,p(m≥p)分別為辛積分傳播子的近似級數(shù)與階數(shù)；cl,dl分別為辛積分傳播子系數(shù)，其取值如表1所示.

中國作為貿(mào)易出口大國，主要特征之一就是出口商品結(jié)構(gòu)低級化。高鐵涉及眾多產(chǎn)業(yè)及技術(shù)，中國高鐵走出國門，就能在極大程度上改善現(xiàn)狀，帶動大量產(chǎn)業(yè)迅速發(fā)展，實現(xiàn)技術(shù)上的不斷突破。另外，龐大的市場必將造就一批擁有高端技術(shù)資金雄厚的國際化企業(yè)，帶動包括高鐵在內(nèi)的一批高科技產(chǎn)品的對外輸出，推動我國出口產(chǎn)品的轉(zhuǎn)型，高端產(chǎn)品將代替低端技術(shù)、低端商品進入國際市場，推進中國出口貿(mào)易由產(chǎn)品貿(mào)易向技術(shù)貿(mào)易的轉(zhuǎn)型，有助于“中國制造”向“中國創(chuàng)造”轉(zhuǎn)變，影響深遠。

表1 辛積分傳播子系數(shù)Tab. 1 The coecients of symplectic integrator propagator

從式(8)可以看出，辛算法具有顯式逐級遞推的特點. 因此，同隱式算法和RK-MRTD算法相比，能夠有效節(jié)省計算資源，提高計算效率.

在空間方向和時間方向分別利用MRTD算法、辛積分傳播子算法離散后，建立了Maxwell方程的離散辛框架，即S-MRTD算法.

為了減少因為電場和磁場幅度差異引起的誤差，經(jīng)過如下歸一化處理：

歸一化后電場和磁場x方向上的迭代公式如下：

式中：

2 數(shù)值穩(wěn)定特性和數(shù)值色散特性分析

2.1 數(shù)值穩(wěn)定性分析

在文獻[10-11]中詳細分析了5級4階的SMRTD算法穩(wěn)定特性，為了簡單起見，在均勻離散的情況，即?x=?y=?z=?時，穩(wěn)定性常數(shù)可以定義為

CFLmax稱為最大穩(wěn)定性常數(shù). 在算法的使用中，穩(wěn)定性常數(shù)的選取不能超過最大穩(wěn)定性常數(shù)，否則會造成程序的不穩(wěn)定，數(shù)值計算結(jié)果出現(xiàn)發(fā)散. 基于D2小波尺度函數(shù)S-MRTD算法的最大穩(wěn)定數(shù)常數(shù)列于表2中，為了能更好地對比S-MRTD算法對于穩(wěn)定性條件的改善，基于D2小波尺度函數(shù)MRTD算法的最大穩(wěn)定數(shù)常數(shù)也列在表2中. 從表2可以看出，相對于5級4階辛算子，優(yōu)化的3級3階辛算子[12-13]對于算法的數(shù)值穩(wěn)定特性有明顯的改善.

表2 MRTD算法的最大穩(wěn)定性常數(shù)Tab. 2 The CFLmax number of MRTD method

2.2 數(shù)值色散特性分析

使用統(tǒng)一的空間步，穩(wěn)定性常數(shù)取CFL=0.4來繪制色散曲線. 圖1(a)為算法相對相速度誤差隨空間分辨率(points per wavelength, PPW)的變化，圖1(b)為算法相對相速度誤差隨穩(wěn)定性常數(shù)的變化，其中平面波的入射角度為 θ=60?，φ=45?.

圖1 相對相速度誤差隨PPW和CFL變化Fig. 1 Relative phase-velocity error vs. PPW and CFL

從圖1(a)可以發(fā)現(xiàn)，3級3階的S-MRTD算法和5級4階的S-MRTD算法都具有良好數(shù)值色散特性，幾乎和RK-MRTD算法相當，明顯比ADI-MRTD算法和MRTD算法更為優(yōu)越，且在相同的計算條件下，S-MRTD算法的數(shù)值計算結(jié)果精度更高. 從圖1(b)可以看到，3級3階的S-MRTD算法和5級4階的S-MRTD算法都具有良好的數(shù)值穩(wěn)定特性，同時有比RK-MRTD算法更為寬松的穩(wěn)定性條件，這明顯改善了MRTD算法苛刻的穩(wěn)定性條件限制，而且在較高的穩(wěn)定常數(shù)條件下，S-MRTD算法也能得到精確的數(shù)值計算結(jié)果. 由于3階3級的SMRTD算法相較于5級4階的S-MRTD算法有更少的時間步，更高的計算效率，故后續(xù)算例皆采用3階3級的S-MRTD算法.

3 數(shù)值計算結(jié)果與分析

3.1 填充介質(zhì)的諧振頻率

矩形波導諧振腔的結(jié)構(gòu)圖同參考文獻[1]，尺寸a×b×c=1m×2m×1.5m ，空間步長 ?s=0.25 m，穩(wěn)定度性數(shù)CFL=0.4，時間步長 ?t=3.3×10?10s，總計算時間步長為50000步，調(diào)制脈沖源頻譜范圍為15～40 MHz，掃描中心頻率f0=30 MHz，掃頻間隔0.01 MHz. 矩形波導諧振腔的頻譜分布如圖2所示. 諧振頻率的解析解與仿真計算結(jié)果見表3，可見S-MRTD算法具有良好的計算精度.

圖2 矩形波導諧振腔的頻譜分布Fig. 2 Spectrum distribution of rectangular waveguide cavity

表3 諧振頻率解析解與仿真結(jié)果對比Tab. 3 Comparison between analytical solution of resonance frequency and simulation results

3.2 波導的截止頻率

非填充介質(zhì)WR-3矩形波導，波導的橫向尺寸a×b=0.864mm×0.432 mm，空間步長?s=0.072 mm，穩(wěn)定性常數(shù)CFL=0.4, 時間步長 ?t=0.096 ps，總計算時間步長為10000步，調(diào)制脈沖源頻譜范圍為170～800 GHz，掃描中心頻率f0=475 GHz，掃頻間隔為0.5 GHz. 在頻譜范圍內(nèi)激勵的模式為 T E10、T E01、TE20、TE11(TM11)、TE21(TM21)、TE12(TM12)、TE22(TM22)，各激勵模式的截止頻率如圖3所示. 歸一化電場能量其峰值對應(yīng)的頻率就是矩形波導的截止頻率，從圖3可以看出，S-MRTD算法對應(yīng)的數(shù)值截止頻率與解析解吻合較好. 在相同誤差邊界(0.3%)條件下，改變空間步長和CFL，MRTD算法和S-MRTD算法仿真所消耗的內(nèi)存和CPU時間見表4. 從表4中數(shù)據(jù)可以看出，S-MRTD算法比MRTD算法大約節(jié)省69.1%的內(nèi)存和58.5%的計算時間.

圖3 WR-3矩形波導頻譜范圍內(nèi)各激勵模式的截止頻率Fig. 3 The cut off frequency of rectangular waveguide

表4 兩種算法消耗的計算機資源對比Tab. 4 The consumption of 2 computation resources

3.3 填充介質(zhì)的波導散射參數(shù)

部分填充介質(zhì)的WR-90矩形波導結(jié)構(gòu)如圖4所示, 波導橫向尺寸a×b=22.86mm×10.16 mm，介質(zhì)塊放置在矩形波導底面上，空間步長?s=1.27 mm，時間步長 ?t=1.693 ps，穩(wěn)定性常數(shù)CFL=0.4，總計算時間步長為20000步，調(diào)制脈沖源頻譜范圍為8～13 GHz，掃描中心頻率f0=10 GHz，掃頻間隔0.03125 GHz. 為了使S-MRTD算法在粗網(wǎng)格和較高的穩(wěn)定性常數(shù)下仍能得到精確的計算結(jié)果，參數(shù)設(shè)置如下：空間步長 ?s=2.54 mm，穩(wěn)定性常數(shù)CFL=0.75，圖5為矩形波導散射參數(shù)的數(shù)值結(jié)果.從圖5可以看出，本文計算的散射參數(shù)與高階辛時域有限差分(symplectic finite-difference time-domain，S-FDTD)算法[14]計算結(jié)果吻合度很好.

圖4 內(nèi)部填充介質(zhì)的WR-90矩形波導Fig. 4 WR-90 rectangular waveguide filled with dielectric material

圖5 WR-90矩形波導散射參數(shù)結(jié)果對比Fig. 5 Comparison of scattering parameters of WR-90rectangular waveguide

4 結(jié) 論

本文把優(yōu)化的3級3階辛算子引入到S-MRTD算法，對其數(shù)值穩(wěn)定性和數(shù)值色散特性的分析顯示了其對于ADI-MRTD和RK-MRTD算法在數(shù)值穩(wěn)定和數(shù)值色散方面的優(yōu)化特性. 應(yīng)用該算法對矩形波導結(jié)構(gòu)進行了研究，相關(guān)特征參數(shù)的數(shù)值仿真結(jié)果顯示了優(yōu)化S-MRTD算法的正確性，并且消耗更少的CPU時間和計算機內(nèi)存. S-MRTD算法的引入為波導相關(guān)特性的分析提供了新的高效的計算算法.