系統(tǒng)矩陣表示的電磁場時域計算方法

鄭宏興 張玉賢

（天津職業(yè)技術(shù)師范大學 天線與微波技術(shù)研究所，天津300222）

引 言

麥克斯韋電磁學理論體系的建立，對現(xiàn)代文明產(chǎn)生了重大影響.伴隨計算機技術(shù)的快速發(fā)展，人們采用有限元法［1-2］、矩量法［3-4］和時域有限差分（Finite Difference Time Domain，F(xiàn)DTD）法［5-6］等數(shù)值計算方法，通過編程來求解麥克斯韋方程組，模擬空間中的電磁場分布，從而解決復雜結(jié)構(gòu)的電磁學工程應(yīng)用問題，目前已經(jīng)廣泛應(yīng)用于微波與毫米波通信、雷達、導航和地質(zhì)探測等各個領(lǐng)域.FDTD是在時域內(nèi)直接求解麥克斯韋方程組，標志性的工作是Yee提出的立方體網(wǎng)格剖分［5］，電場和磁場分開半個時間步交錯迭代，優(yōu)點是計算的物理過程直觀，對復雜目標建模相對容易，在計算機內(nèi)部存儲變量按照時間遞推過程進行更新，對計算機內(nèi)存要求相對較低，但是時間步長和空間網(wǎng)格受到穩(wěn)定性條件的限制［6］.后 來，Namiki［7］和Zheng［8］分 別 提 出了無條件穩(wěn)定技術(shù)，從而突破了時間離散間隔的限制，采用變換方向的隱式（Alternative Direction Implicit，ADI）方法來求解方程，可實現(xiàn)無條件穩(wěn)定.為了動態(tài)表達上述數(shù)值計算結(jié)果，一方面需要應(yīng)用計算機的圖形顯示功能直觀顯示，另一方面還需要應(yīng)用數(shù)字信號處理（Digital Signal Processing，DSP）［9-12］技術(shù)表達電磁波與物質(zhì)的相互作用過程.DSP也是一門伴隨計算機技術(shù)發(fā)展的數(shù)值計算技術(shù)，在語音、圖像、地質(zhì)探測和生物醫(yī)學等方面得到廣泛應(yīng)用.無論是把DSP與計算電磁學相結(jié)合，用來處理等離子體或色散媒質(zhì)的時間移位算子問題［9-10］，還是將N階色散媒質(zhì)的研究與DSP技術(shù)相結(jié)合求解復雜結(jié)構(gòu)問題［11］，或者用系統(tǒng)函數(shù)來研究時域電磁問題等［12］，這些都說明了兩種方法相結(jié)合對電磁場數(shù)值計算的深入研究是很有意義的.

為了充分利用DSP技術(shù)的研究成果，本文按照信號與系統(tǒng)理論的概念，把電磁場的求解區(qū)域看做一個線性系統(tǒng)，討論關(guān)于這個系統(tǒng)方程的建立過程.將電磁場方程的時域離散差分方程［13-14］轉(zhuǎn)換為系統(tǒng)矩陣方程，通過分析時域迭代過程，將其納入系統(tǒng)仿真問題［13-15］的框架下，結(jié)合DSP方法來實現(xiàn)將多維空間的多變量運算轉(zhuǎn)換為一維矢量的運算，進一步給出在無條件穩(wěn)定下電磁場方程所對應(yīng)的系統(tǒng)矩陣的表示形式.利用系統(tǒng)矩陣求解麥克斯韋方程組，對時域和空域的電磁環(huán)境進行模擬仿真研究，達到分析電磁場問題、提高計算效率的目的.

1 電磁場方程的時域離散矩陣形式

考慮媒質(zhì)的一般特性，有源區(qū)域中的麥克斯韋方程組可以表達為如下的對稱形式［16］：

式中：Je＝σe·E，Jm＝σm·H分別為電流源和等效磁流源，σe和σm分別為電導率和等效磁損耗張量.以先計算磁場H，后計算電場E為主要計算過程，故設(shè)磁場H的時間步在n＋1／2，電場E的時間步在n＋1，如此實現(xiàn)式（1）和（2）的時域離散差分格式，得到：

式中，C為旋度算子▽×的矩陣形式，定義為

它是一個反對稱張量.式（3）和（4）為麥克斯韋方程組的時域離散矩陣形式，按照時間遞推，我們可以得到空間上波的傳播過程.

2 時域差分方法的系統(tǒng)方程推導

由式（3）和（4），容易得到時域離散的電磁計算過程所對應(yīng)的系統(tǒng)框圖，如圖1所示.

圖1 傳統(tǒng)FDTD計算過程的系統(tǒng)框圖

顯然，時域離散的電磁計算方法是利用時間上的交錯迭代性質(zhì)來模擬整個空間上的電磁場分布，并以此來描述空間中的電磁波傳播過程.從式（4）發(fā)現(xiàn)，En＋1中存在由Hn＋1／2決定的場分量.將式（4）展開，

將式（3）中的磁場分量H移動1／2個時間步，得

設(shè)x［n］＝｛［E H］n｝T是由電場E和磁場H在n時刻所有空間離散后的各點場量.如果把求解區(qū)域看做線性系統(tǒng)，那么式（6）和（7）描述的系統(tǒng)在n＋1時刻的電磁場分布x［n＋1］，可以由圖2所示的系統(tǒng)框圖表達，圖中R描述的是線性系統(tǒng)的系統(tǒng)矩陣，定義為

圖2 關(guān)于系統(tǒng)矩陣R的電磁場迭代原理圖

由式（6）和（7）分別得到：

由系統(tǒng)矩陣R中的第一項子矩陣K可見，重復對電場E進行了兩次空間微分運算，C矩陣的元素為空域上一階導數(shù)，作用兩次后使得電場E在空域上進行了二階求導.從這里開始，在式（3）和（4）中所考慮的是多變量的交替問題，均可歸結(jié)到數(shù)字信號處理的問題，其優(yōu)勢在于多變量的問題直接轉(zhuǎn)換為一維矢量的問題，重點在于系統(tǒng)矩陣R中的元素分布情況.考慮到差分方程的求解過程中，這種電磁計算的系統(tǒng)方程應(yīng)滿足：

式中：S為激勵源分布矩陣；u［n］為輸入的激勵信號；y［n］為輸出信號；c，d均為常數(shù).從式（3）和（4）轉(zhuǎn)換為式（10）和（11），場方程的時域矩陣形式直接變換為一維矢量的運算，所有的電磁場量由一維矢量直接反映出來，該方法在不同維數(shù)下的存儲變量數(shù)與傳統(tǒng)FDTD方法的對比結(jié)果如表1所示.在滿足穩(wěn)定性條件的前提下，系統(tǒng)矩陣方法比傳統(tǒng)FDTD的存儲變量數(shù)目明顯減少，能有效提高計算機的效率.

表1 不同維數(shù)下系統(tǒng)矩陣方法與傳統(tǒng)FDTD的存儲變量數(shù)

3 無條件穩(wěn)定的系統(tǒng)矩陣形式

對于式（3）和（4），依據(jù)時域離散差分格式來進行處理.由于空間微分的矩陣算子C所對應(yīng)的場量E和H均為同時刻，導出式（10）和（11）的結(jié)果與傳統(tǒng)FDTD一樣，都需要滿足時間離散間隔的穩(wěn)定性條件.對于無條件穩(wěn)定的FDTD［7-8］，一個計算時間步分為兩個過程進行，這里對式（5）中的微分算子C所作用的場量，分別進行時域上的交替處理，需要兩個計算過程來描述.為了便于矩陣運算的形式表達，將C分解為

式中：

在過程一中，磁場H的時間步在n＋1／4，電場E的時間步在n＋1／2；過程二中，磁場H的時間步在n＋3／4，電場E的時間步在n＋1.可得到：

過程一

過程二

在 過 程 一 中，Hn＋1／4和En＋1／2都 是 未 知 的，Hn-1／4和En都是已知的，不妨將式（14）代入式（15）中，并將En＋1／2的所有項放在等式的左邊，令

其中I為單位矩陣，可得到

同樣，將式（16）代入式（17）中，導出En＋1的關(guān)系式：

根據(jù)式（14）、（18）和（16）、（19）可分別得到對應(yīng)兩個計算過程的系統(tǒng)框圖，如圖3所示.

圖3 無條件穩(wěn)定FDTD系統(tǒng)框圖

無條件穩(wěn)定的FDTD系統(tǒng)方程組表示相對較復雜，不妨假設(shè)：

式中K1，L1和K2，L2分別為兩個過程對應(yīng)的系統(tǒng)矩陣R的第一行子矩陣，其表示為：

于是式（14）和（16）變換為：

分別將式（23）和（24）中的未知量展開，得：

式中M1，N1和M2，N2分別為兩個過程對應(yīng)的系統(tǒng)矩陣R的第二行子矩陣，其表示為：

因此，無條件穩(wěn)定的FDTD方程組對應(yīng)的系統(tǒng)矩陣的元素可由式（20）、（21）和（25）、（26）分別表示.由于在各向同性媒質(zhì)中，Φ1和Φ2均為三條帶矩陣，它的求逆運算方法在計算機技術(shù)中相當成熟，采用系統(tǒng)矩陣表示的形式，非常便于計算機的編程.

4 數(shù)值實驗

4.1 一維高斯脈沖的傳播

將式（3）和（4）的結(jié)果投影在一維空間上進行離散化處理，可得到：

式中

將上述參數(shù)初始化，在一維計算區(qū)域的中心，設(shè)置高斯脈沖形式的電流源為

式中：t0＝2×10-10；τ＝5×10-11.初始化設(shè)置Δx＝1mm，如圖4所示，圖中x＝0m和x＝1m平面為導電壁，實線為電場，虛線為磁場，從圖中觀察電磁波在一維系統(tǒng)矩陣方法下的仿真過程.電流源從中心處向左右兩邊傳播，碰到導體面后產(chǎn)生反射.

圖4 不同時間步的系統(tǒng)矩陣方法動態(tài)過程

為了對比系統(tǒng)矩陣方法與傳統(tǒng)FDTD的計算效率，采用主頻為2.8GHz、雙核CPU、內(nèi)存為4G比特的計算機進行仿真，迭代次數(shù)分別設(shè)置為1 000和5 000，其計算時間如表2所示.可知，相對傳統(tǒng)FDTD，系統(tǒng)矩陣方法的計算效率提高23%.

表2 兩種電磁計算方法的時間對比

4.2 低通濾波器

通過式（20）、（21）以及式（25）、（26），得到了在無條件穩(wěn)定的系統(tǒng)矩陣方法.s為CFL（Courant-Friedrichs-Lewy）因子［7］.傳統(tǒng)FDTD方法必須滿足s≤1的條件，但在無條件穩(wěn)定的FDTD計算中，s的取值可以大于1.這里以低通濾波器為例，設(shè)空間網(wǎng)格尺寸分別為Δx＝0.5mm，Δy＝0.5mm，Δz＝0.3mm，選取接地板尺寸為50Δx×46Δy，介質(zhì)板的相對介電常數(shù)為εr＝2.2，其厚度為3Δz＝0.9 mm，結(jié)構(gòu)及坐標網(wǎng)格如圖5所示.應(yīng)用傳統(tǒng)FDTD方法、ADI-FDTD方法以及系統(tǒng)矩陣方法分別對這種結(jié)構(gòu)的低通濾波器進行計算，結(jié)果如表3所示.

圖5 低通濾波器的結(jié)構(gòu)及其尺寸坐標圖

傳統(tǒng)FDTD方法受Courant穩(wěn)定性條件的限制，即使CFL因子s＝1，其時間步仍需要2 700步才能完成計算.ADI-FDTD方法雖然能夠突破Courant穩(wěn)定性條件的限制，但是計算時間反倒增加了，這種方法并不能減少仿真時間.本文的系統(tǒng)矩陣方法，不受Courant穩(wěn)定性條件的限制，當s＝3時，與傳統(tǒng)FDTD方法相比，減少了一半以上的計算時間.圖6給出三種方法的仿真結(jié)果，系統(tǒng)矩陣方法與ADI-FDTD的結(jié)果吻合.同時發(fā)現(xiàn)當s＝0.9時，它們的結(jié)果幾乎重合.因ADI-FDTD方法隨CFL因子s的增大，Δt也就增大，對應(yīng)的頻率就會減小，因此其頻譜響應(yīng)將會向低頻方向移動，圖6的散射參數(shù)S11曲線也會向低頻方向移動.說明本文的系統(tǒng)矩陣方法是有效的.

表3 三種數(shù)值方法計算低通濾波器的時間

圖6 低通濾波器的S11參數(shù)

4.3 緊湊型帶阻濾波器

為了進一步驗證系統(tǒng)矩陣方法的有效性，選取D.Cheng在文獻［17］中提出的緊湊型帶阻濾波器，其結(jié)構(gòu)如圖7所示，其中W3＝0.7mm，G＝W4＝W5＝0.5mm，L3＝24.2mm，L4＝21.1mm和L5＝20.6mm.這種結(jié)構(gòu)對于x方向的網(wǎng)格剖分需要更加精細.考慮到L4和L5的長度相近，構(gòu)成的縫隙至少需要兩個網(wǎng)格以上，于是設(shè)置網(wǎng)格尺寸為Δx＝0.54mm，Δy＝0.768mm，Δz＝0.114mm，分別采用傳統(tǒng)FDTD和本文的方法進行計算.為了更好地說明它的實用性，用商業(yè)軟件HFSS15.0.2仿真，并將其結(jié)果與上述兩種方法對比，如圖8所示.可知，在CFL因子s取值稍大時，緊湊型帶阻濾波器的散射參數(shù)S11頻率略有移動，但計算結(jié)果依然滿足要求.說明本文方法與商業(yè)軟件HFSS15.0.2的仿真結(jié)果相媲美.

圖7 緊湊型帶阻濾波器的結(jié)構(gòu)圖

圖8 緊湊型帶阻濾波器的S11參數(shù)

5 結(jié) 論

采用系統(tǒng)矩陣方法來研究時域電磁問題，是充分應(yīng)用DSP的研究成果所提供的計算機軟件資源的一次嘗試.對于復雜媒質(zhì)的電磁結(jié)構(gòu)，如復雜目標散射、電磁兼容設(shè)計等，針對具體的物理模型，利用DSP技術(shù)建立系統(tǒng)矩陣是一個相對容易的工作.通過系統(tǒng)傳遞函數(shù)的拉普拉斯反變換，建立時域系統(tǒng)矩陣，從而求解這個系統(tǒng)方程，是今后進一步開展的工作.上述算例均表明本文的方法是有效的，能夠顯著提高計算效率.

［1］吳 霞，周樂柱.時域有限元在計算電磁問題上的應(yīng)用及發(fā)展［J］.電波科學學報，2008，23（6）：1209-1216.WU Xia，ZHOU Lezhu.Application and development of time-domain finite element method on EM analysis［J］.Chinese Journal of Radio Science，2008，23（6）：1209-1216.（in Chinese）

［2］JIN J M，LOU Z，LI Y J，et al.Finite element analysis of complex antennas and arrays［J］.IEEE Transactions on Antennas and Propagation，2008，56（8）：2222-2240.

［3］唐 璞，王 建，王超然，等.曲線分段的矩量法分析阿基米德平面螺旋天線［J］.電波科學學報，2006，21（6）：929-934.TANG Pu，WANG Jian，WANG Chaoran，et al.Analysis of archimedean spiral antenna using moment method with curved segmentations［J］.Chinese Jour-nal of Radio Science，2006，21（6）：929-934.（in Chinese）

［4］張小林，高火濤.近地短波魚骨天線矩量法建模及應(yīng)用［J］.電波科學學報，2007，22（2）：281-285.ZHANG Xiaolin，GAO Huotao.Moment method modeling and application for shortwave fishbone antennas above ground［J］.Chinese Journal of Radio Science，2007，22（2）：281-285.（in Chinese）

［5］YEE K S，Numerical solution of initial boundary value problem involving Maxwell’s equations in isotropic media［J］.IEEE Transactions on Antennas and Propagation，1966，14（3）：302-307.

［6］葛德彪，閆玉波.電磁波時域有限差分方法［M］.3版.西安：西安電子科技大學出版社，2011.

［7］NAMIKI T.A new FDTD algorithm based on alternating direction implicit method［J］.IEEE Transactions on Microwave Theory and Techniques，1999，47（10）：2003-2007.

［8］ZHENG F，CHEN Z，ZHANG J.A finite-different time-domain method without the courant stability condition［J］.IEEE Microwave Guided Wave Letter，1999，9（11）：441-443.

［9］張玉強，葛德彪.一種基于數(shù)字信號處理技術(shù)的改進通用色散介質(zhì)移位算子時域有限差分方法［J］.物理學報，2009，58（12）：8243-8248.ZHANG Yuqiang，GE Debiao.An improved shift operator finite-difference time-domain method based on digital signal processing technique for general dispersive medium［J］.Acta Physica Sinica，2009，58（12）：8243-8248.（in Chinese）

［10］葛德彪，吳躍麗，朱湘琴.等離子體散射FDTD分析的移位算子方法［J］.電波科學學報，2003，18（4）：359-362.GE Debiao，WEI Yueli，ZHU Xiangqin.Shift operator method applied for dispersive medium in FDTD analysis［J］.Chinese Journal of Radio Science，2003，18（4）：359-362.（in Chinese）

［11］汪 彤，張文俊，劉維亮.N階色散媒質(zhì)的FDTD法與數(shù)字信號處理技術(shù)［J］.電波科學學報，2000，15（1）：455-458.WANG Tong，ZHANG Wenjun，LIU Weiliang.FDTD and digital signal processing for the Nth-order dispersive media［J］.Chinese Journal of Radio Science，2000，15（1）：455-458.（in Chinese）

［12］陳春雨，吳 群.基于系統(tǒng)函數(shù)的FDTD特性研究［J］.微波學報，2010，26（5）：53-55.CHEN Chunyu，WU Qun.Study on characteristics of FDTD based on system function［J］.Journal of Microwaves，2010，26（5）：53-55.（in Chinese）

［13］ALAN V O，ALAN S W，HAMID N，et al.Signals and Systems［M］.2nd ed.New York：Pearson Education，2013.

［14］吳迎年，張 霖，張利芳，等.電磁環(huán)境仿真與可視化研究綜述［J］.系統(tǒng)仿真學報，2009，21（20）：6332-6338.WU Yingnian，ZHANG Lin，ZHANG Lifang，et al.Survey on electromagnetic environment simulation and visualization［J］.Journal of System Simulation，2009，21（20）：6332-6338.（in Chinese）

［15］殷 勤，陳 彬，汪 瑩，等，電磁跟蹤系統(tǒng)磁傳感器性能的數(shù)值模擬研究［J］.傳感技術(shù)學報，2010，23（8）：1079-1083.（in Chinese）YIN Qin，CHEN Bin，WANG Ying，et al.Study on numerical simulation of electromagnetic tracker system magnetic sensor［J］.Chinese Journal of Sensors and Actuators，2010，23（8）：1079-1083.（in Chinese）

［16］王一平.工程電動力學（修訂版）［M］.西安：西安電子科技大學出版社，2005.

［17］CHENG D，YIN H C，ZHENG H X.A compact dual-band bandstop filter with defected microstrip slot［J］.Journal of Electromagnetic Waves and Applications，2012，26（10）：1374-1380.