用于各向異性媒質(zhì)的TLM算法

夏新仁，張國華，金賢龍，謝 輝

（空軍駐合肥地區(qū)軍事代表室，安徽 合肥 230000）

0 引言

近年來，各向異性媒質(zhì)的電磁問題受到關(guān)注，出現(xiàn)大量用于處理各向異性媒質(zhì)的算法，特別是用于處理各向異性媒質(zhì)的改進FDTD算法。1990年，LUEBBERS提出一種用離散卷積處理色散媒質(zhì)的RC－FDTD算法，并將此算法推廣到各向異性的磁化等離子體，處理各向異性的磁化鐵氧體及手征媒質(zhì)［1－3］。1996年，KELLEY采用電場的分段線性近似——PLRC－FDTD改善了RC－FDTD算法的計算精度［4］。處理各向異性媒質(zhì)的FDTD方法還有JEC法、DI法、分段線性電流密度遞歸卷積（PLJERC）法、ADE法等［5－8］。由于整個網(wǎng)格區(qū)間的離散化使所需的網(wǎng)格數(shù)巨大，盡管FDTD的計算復(fù)雜度相對較低，但所需的存儲時間和CPU時間仍較長。

TLM算法是一種用于處理電磁問題的全新時域數(shù)值方法，該方法功能強大且計算精度高。2002年，YAICH運用TLM算法結(jié)合HSCN技術(shù)處理磁化等離子體的電磁問題，并擴展至磁化鐵氧體［9－10］。本文在Z變換中引入TLM算法，對各向異性媒質(zhì)的電磁性質(zhì)計算進行了研究，并與解析法結(jié)果進行了對比［11］。

1 各向異性媒質(zhì)的TLM算法

對一般各向異性媒質(zhì)，麥克斯韋方程可寫成

設(shè)電流密度和磁流密度為零，一般各向異性媒質(zhì)的本構(gòu)關(guān)系為

式中：H為磁場強度；E為電場強度；ε0為自由空間介電常數(shù)；μ0為真空磁導(dǎo)率；“＊”為時域卷積；，分別為電極化率和磁化率張量，兩者均為3×3維張量的時間函數(shù)。本文研究一般各向異性媒質(zhì)，非雙各向異性媒質(zhì)，故式（2）中的耦合項為零。

將式（2）代入式（1），可得

假設(shè)電磁波沿x向傳播，以一維電磁問題為例，在電磁波傳播的x方向無場分量，y、z向存在相互耦合的電磁場分量，場分量空間節(jié)點如圖1所示［12］。

圖1 一維各向異性媒質(zhì)的TLM節(jié)點Fig.1 TLM node of one－dimensional isomerism medium

場分量空間節(jié)點為一個邊長為Δl的立方體，4個場分量Ey，Ez，Hy，Hz設(shè)置在節(jié)點中心，連接相鄰節(jié)點的端口為V4，V5，V10，V11，端口編號與三維空間一致。因此，在一維空間，電磁波沿x向傳播，式（3）可轉(zhuǎn)化為

式中：s為?／?t轉(zhuǎn)變Laplace域生成；（▽×H）y為▽×H的y方向分量，其余旋度運算下標符號意義相同。電、磁極化率張量可表示為

將式（4）用場等效電路參數(shù)表示，可得

將式（6）變換成行波的形式為

此處，上標“i”為入射波量［11］。定義式（7）左邊為外部激勵，則

定義傳輸矩陣t＝2（2I＋sΔtM）－1，則式（9）轉(zhuǎn)化為

對式（9），用雙線性Z變換可得離散時域，其中sΔt對應(yīng)到Z域為

由部分分式展開可得

式中：M2為z的函數(shù)；M0，M1為常數(shù)張量。

將式（12）代入式（9），整理可得

式中：T－1＝2I＋2M0；κ＝－2I－2M1。

式（13）可簡化為

式中：S＝2Fe＋κF－2M2（z）F。

Z域中的z－1在時域?qū)?yīng)為一個時延算子，式（14）轉(zhuǎn)變到時域可得

對一般媒質(zhì)，傳輸矩陣t因媒質(zhì)的頻率特性不同而改變。用迭代過程計算反射端口電壓，即下一個時間步長鄰近節(jié)點的端口電壓，即

對一維各向異性媒質(zhì)的TLM節(jié)點，傳播速度

式中：Δl為TLM空間節(jié)點的邊長；Δt為時間步長。

對磁化等離子體，若電磁波沿x向傳播，且假設(shè)背景磁場方向平行于x向，則磁化等離子體的等效相對介電常數(shù)張量可表示為

式中：ε1＝1－；ε2＝－；ε3＝1－；j為虛數(shù)單位。此處：νc為電子的碰撞頻率；ω為電磁波頻率；ωp為等離子體頻率；ωb為電子回旋頻率，且ωb＝eB0／me。其中：B0為背景磁場強度；e，me分別為電子電荷數(shù)和質(zhì)量［12］。

電磁波沿x向傳播，電磁場在x向無耦合，則

式 中：χ1＝－；χ2＝－j。

為運用TLM算法，將χ1、χ2由頻域變換到Laplace域，由Laplace域可變換到Z域。將χ1變換到Laplace域可得

運用部分分式展開，則式（20）可轉(zhuǎn)化為

式（21）變形可得

同理，將χ2變換到Laplace域可得

1／s由Laplace域變換到Z域為

用1－z－1乘式（24）右部可得

為得到式（12）形式，用1－z－1乘式（27）的右部，并用部分分式展開可得

將θ0、θ1、θ2乘以－j可將變換到Z域。用θ0，θ1，θ2，α0，α1可構(gòu)建矩陣M0，M1，M2（z）。因此，可得（1－z－1）χ1（s）的Z域形式為

式中：n＝θ0＋z－1（θ1＋θ2）。

式中：m＝－θ1－jB1θ1。

對式（30）進行替代、簡化可得

將（1－z－1）χ2（s）變換到Z域形式為

對式（32）進行替代、簡化可得

構(gòu)建矩陣

由式（13）、（14）可知

令M2（z）F＝Se，由M2（z），F(xiàn)可知Se的第一個元素，即Re［Se1］為

將θ2＝代入式（38），整理可得

由式（39）可得時域的迭代關(guān)系為

同理可得Se的第二個元素，即?［Se2］的時域的迭代關(guān)系為

因此，

2 仿真

為檢驗TLM算法在各向異性媒質(zhì)中的有效性和計算精度，對某一維實例用TLM算法模擬電磁脈沖垂直入射到各向異性磁化等離子體層中的傳播。假設(shè)電磁波的傳播方向和外加磁場的方向平行，等離子體層的兩邊為自由空間，入射源采用高斯脈沖。TLM計算空間分為500個空間節(jié)點元胞，空間節(jié)點編號0～300和420～500為真空，等離子體所占的空間節(jié)點編號取300～420，空間步長取Δ＝75μm，假設(shè)條件穩(wěn)定，時間步長取Δt＝0.125ps。背景磁場強度B0＝1.7T，等離子體密度n0＝3.1×1019m－3，碰撞頻率νc＝20GHz，磁化等離子體層厚度d＝9mm。

電磁波在各向異性磁化等離子體中傳播時，其基本模式為左旋圓極化波（LCP）和右旋圓極化波（RCP）。為得到LCP和RCP的反射系數(shù)與透射系數(shù)的幅度及相位與頻率的變化關(guān)系，實現(xiàn)過程為：對等離子體層的前后的電場強度值作離散Fourier變換［13］。將電場強度隨時間的變化轉(zhuǎn)換到頻域，得到電場隨頻率的變化關(guān)系，有

用解析法和TLM算法計算各向異性磁化等離子體層中左、右旋圓極化波的反射系數(shù)和透射系數(shù)的幅度和相位結(jié)果分別如圖2～9所示［14］。由圖2～9可知：本文結(jié)果與解析的結(jié)果吻合得較好，表明TLM算法具有高效性和高精度。

圖2 右旋極化波反射系數(shù)幅值Fig.2 Reflection coefficient magnitude for RCP polarization

圖3 右旋極化波反射系數(shù)相位Fig.3 Reflection coefficient phase for RCP polarization

圖4 左旋極化波反射系數(shù)幅值Fig.4 Reflection coefficient magnitude for LCP polarization

圖5 左旋極化波反射系數(shù)相位Fig.5 Reflection coefficient phase for LCP polarization

圖6 右旋極化波透射系數(shù)幅值Fig.6 Transmission coefficient magnitude for RCP polarization

圖7 右旋極化波透射系數(shù)相位Fig.7 Transmission coefficient phase for RCP polarization

圖8 左旋極化波透射系數(shù)幅值Fig.8 Transmission coefficient magnitude for LCP polarization

3 結(jié)束語

本文提出了一種用于各向異性媒質(zhì)的全新時域數(shù)值方法。該法可對各向異性媒質(zhì)進行建模，也可用于各向異性媒質(zhì)中波的傳播和散射分析。方法通過Z變換由頻域變換到Z域，避免了時域的卷積運算，保證了計算的高效率。用本文方法計算了磁化等離子體平板對平行于磁場傳播的電磁波的反射和透射系數(shù)幅度和相位，并與解析解的結(jié)果進行比對，兩種方法吻合較好，驗證了方法的高效性和高精度性。

圖9 左旋極化波透射系數(shù)相位Fig.9 Transmission coefficient phase for LCP polarization

［1］LUEBBERS R，HUNSBERGER F P，KUNZ K S，et al.A frequency－dependent finite－difference time－domain formulation for dispersive materials［J］.Transactions on Electromagnetic Compatibility，1990，32（3）：222－227.

［2］HUNSBERGER F，LUEBBERS R，KUNZ K.Finite－difference time－domain analysis of gyrotropic media magnetized plasma［J］.Transactions on Antennas and Propagation，1992，40（12）：1489－1495.

［3］HUNSBERGER F P.Extension of the finite－difference time－domain method to gyrotropic media［D］.Nittany Valley：The Pennsylvania State University，1991.

［4］KELLEY D，LUEBBERS R J.Piecewise linear recursive convolution for dispersive media using FDTD［J］.Transactions on Antennas and Propagation，1996，44（6）：792－797.

［5］XU L J，YUAN N C.JEC－FDTD for 2－D conducting cylinder coated by anisotropic magnetized plasma［J］.IEEE Microwave and Wireless Components Letters，2005，15（12）：892－894.

［6］劉少斌，莫錦軍，袁乃昌.各向異性磁化等離子體的JEC－FDTD算法［J］.物理學(xué)報，2004，53（3）：783－787.

［7］LIU S B，MO J J，YUAN N C.Piecewise linear current density recursive convolution FDTD implementation for anisotropic magnetized plasmas［J］.Microwave and Wireless Components Letters，2004，14（5）：222－224.

［8］劉少斌，莫錦軍，袁乃昌.各向異性磁化等離子體的輔助方程FDTD算法［J］.物理學(xué)報，2004，53（7）：2233－2236.

［9］YAICH M I，KHALLADI M，ZEKIK I，etal.Modeling of frequency－dependent magnetized plasma in hybrid symmetrical condensed TLM method［J］.Microwave and Wireless Components Letters，2002，12（8）：293－295.

［10］YAICH M I，KHALLADI M.A SCN－TLM model for the analysis of ferrite media［J］.Microwave and Wireless Components Letters，2003，13（6）：217－219.

［11］LIU S，ZHONG S Y，LIU S B.Piecewise linear recursive convolution FDTD method for magnetized plasma［J］.Journal of Systems Engineering and Electronics，2006，17（2）：290－295.

［12］FURSE C M，GANDHI O P.Why the DFT is faster than the FFT for FDTD time－to－frequency domain conversions［J］.IEEE Microwave Guided Wave Letter，1995，9（5）：326－328.

［13］袁敬閎，莫懷德.等離子體中的波［M］.成都：電子科技大學(xué)出版社，1990.