有耗介質中時域精細積分方法的數值色散分析

孫 剛 馬西奎 白仲明

（西安交通大學，電氣工程學院，電力設備電氣絕緣國家重點實驗室，陜西 西安710049）

引 言

時域有限差分（FDTD）方法［1－2］是一種有效的電磁場時域計算方法，并已廣泛應用于電磁散射、電磁兼容、波導與諧振腔系統、天線輻射特性的研究［3－6］。然而，FDTD的空間步長要受到數值色散條件的限制，其時間步長要受到柯朗（Courant）穩定性條件的限制。一般而言，空間步長要小于電磁波波長的1／10［2，7］。這些限制條件使得應用 FDTD分析電大尺寸（仿真模型的幾何尺寸遠大于電磁波的波長）問題時所需求的計算機CPU時間和內存過大，出現計算困難。

為了克服FDTD的上述固有缺陷，文獻［8］提出了交替方向隱式時域有限差分（ADI－FDTD）方法。ADI－FDTD是無條件穩定的，但是數值色散誤差較大。文獻［9］的研究表明：當電導率較大時，ADI－FDTD的應用要受到限制。

為了解決FDTD和ADI－FDTD的上述問題，時域精細積分（PITD）方法［10］將精細積分技術［11］引入到時域電磁場計算之中。PITD的時間步長可以遠大于Courant條件所限制的值［12］。在無耗介質中，PITD的數值色散誤差略大于FDTD的數值色散誤差，小于ADI－FDTD的數值色散誤差，且基本上不受時間步長的影響［13］。

文獻［13］對PITD的數值色散分析僅限于無耗介質。分析有耗介質中PITD的數值損耗和數值色散，將為其在有耗介質中的應用和進一步發展提供理論基礎。

1.數值色散方程

1.1 時域精細積分方法的基本原理

在PITD中，首先應用Yee網格對兩個麥克斯韋（Maxwell）旋度方程

進行空間離散，得到如下的常微分方程組［10］：

式中：X＝ ［Ex，Ey，Ez，Hx，Hy，Hz］T是由網格結點上的電場分量和磁場分量構成的列向量，M是一個由空間步長（Δx，Δy，Δz）決定的系數矩陣。若引入空間位移算子Xh，Yh和Zh［14］，就可以將 M 矩陣表示成如下的形式

Dx＝和是空間差分算子。

根據常微分方程組理論，可以得到X的時間遞進公式

式中：Δt是時間步長；T＝exp（Δt M）.

指數矩陣T的高效率、高精度計算是式（4）遞進計算的關鍵。在PITD中，采用如下的有限項泰勒級數展開對T進行精細積分計算［11］

式中：τ＝Δt／2N是子時間步長，N是預先選定的整數，一般取N≥20.有關應用精細積分技術計算式（5）的詳細過程見文獻［10］［11］。

1.2 一維數值色散方程

在一維情況下（X＝ ［Ex，Hy］T），矩陣M 就簡化成

對于沿著z方向傳播的正弦均勻平面波，X在時間－空間譜域內，可以表示成如下的離散形式［2］

其中：γ是數值傳播常數。此時有

不難看出：矩陣M的特征值λM要滿足如下的方程

注意到矩陣M是可對角化的

式中，Y是由矩陣M的特征向量構成的矩陣。將式（10）代入式（5），有

即矩陣T也是可對角化的，并且有

式中，λT是矩陣T的特征值。將式（7）代入式（4），就可以得到

上述方程的非平凡解要求

利用式（11）和式（12），可以將式（14）簡化成

在 PITD 中，τ 是 一 個 很 小 的 量［10－11］，因 此，式（15）的右邊可以近似表示成

比較式（15）和式（16），不難看出：λM≌jω.不妨設λM＝jη，此時式（15）可以表示成

將γ＝α＋jβ（α是數值衰減常數，β是數值相位常數）代入式（9）就可以得到

式（17）和式（18）共同構成了一維情況下有耗介質中PITD的數值色散方程。

1.3 三維數值色散方程

在三維情況下，PITD的數值色散方程可以由一維情況直接拓展得到，不再贅述其推導過程，只給出其結果如下

式中：w＝x，y，z.并且有

在式（19）中：α和β是通過三角函數和雙曲函數耦合的，這些函數的變量中還包括了電磁波的傳播方向。因此，無法直接計算α和β.α和β的值，但可以應用Newton－Rahphson方法迭代得到。

為了研究數值損耗和數值色散的各向異性，通常要考慮數值波沿坐標軸方向（可以等效為一維情況）和數值網格對角線方向的傳播特性。在一維情況、二維正方形網格情況（Δx＝Δy＝Δ）和三維正立方體網格情況（Δx＝Δy＝Δz＝Δ）下，式（19）可以簡化如下

式中：tanδ＝（ω／η）tanδ0是數值損耗角正切，tanδ0是損耗角正切表示空間維數。由η≌ω可知：tanδ≌tanδ0，即PITD的數值損耗角正切是損耗角正切的高精度近似。式（21）中的α和β是可以進一步分離的，有

式（22）可以用于直接求解數值波沿坐標軸方向和對角線方向傳播情況下的α和β.

對比式（19）和有耗介質中FDTD的數值色散方程式（23）［15］，

可以看出：當時間步長趨近于0時，式（19）和式（23）將趨于一致。這是因為這兩種方法在空間上都采用了Yee網格離散。當時間步長趨近于0時，數值色散方程體現出的就是空間離散網格的數值色散特性。

當電導率σ趨近于0時，式（19）將變為無耗介質中PITD的數值色散方程［13］

當空間步長和時間步長均趨近于0時，式（19）將趨近于有耗介質中電磁波傳播的色散方程

式中：α0是電磁波的真實衰減常數；β0是電磁波的真實波數。

2.數值損耗和數值色散分析

首先分析在一維情況下（電磁波沿x方向傳播，ε＝ε0，μ＝μ0），時間步長 Δt、空間步長 Δx 和電導率σ對PITD的數值損耗和數值色散的影響。然后分析在高維情況下，PITD的數值損耗和數值色散的各向異性。定義數值損耗誤差（NLE）和數值相位誤差（NPE）分別如下

2.1 時間步長的影響

在一維情況下，PITD的NLE和NPE隨時間步長的變化如圖1所示。電磁波頻率f＝300 MHz，空間采樣密度（電磁波波長／空間步長）Nλ＝20，Courant常數S＝cΔt／Δx∈（0，1024），預先選定的整數N＝20，電導率σ分別等于3E－4S／m，3E－2S／m和30S／m（通常電介質材料的電導率σ∈（1E－4，30）S／m）.

在圖1中，所有的曲線都幾乎是平直的。這說明PITD的數值損耗和數值色散都基本上不受時間步長的影響。其原因是PITD的數值色散方程式（19）中所包含的時間量（子時間步長τ）被限制在很小的范圍內（τ≤2.722E－13s）.

2.2 空間步長的影響

在一維情況下，PITD的NLE和NPE隨時間步長的變化如圖2所示。電磁波頻率f＝300 MHz，Courant常數S＝1 024，預先選定的整數N＝20，空間采樣密度Nλ∈（10，40），電導率σ分別等于3E－4S／m，3E－2S／m和30S／m.

在圖2中，隨著空間采樣密度的增大（空間步長減小），NLE和NPE都將趨近于0.這說明空間步長越小，PITD的數值損耗和數值色散的誤差越小。

2.3 電導率的影響

在一維情況下，PITD和FDTD的NLE和NPE隨電導率σ的變化如圖3所示。電磁波頻率f＝300MHz，在PITD中Courant常數S＝1 024，預先選定的整數N＝20.在FDTD中Courant常數S＝1.空間采樣密度Nλ＝20，電導率σ∈（1E－4，30）S／m.

在圖3（a）中，PITD的NLE大于0，即PITD的數值損耗大于電磁波的真實損耗。隨著電導率σ增大，PITD的NLE逐步減小，并和FDTD的NLE趨近于同一極限值。

在圖3（b）中，當電導率較小時（σ?ωε），PITD的NPE大于0.當電導率較大時（σ≥ωε），PITD的NPE小于0，即PITD的數值波速大于電磁波的真實波速。此時，PITD的NPE的絕對值要小于FDTD的NPE的絕對值，即PITD的數值波速比FDTD的數值波速更加接近于電磁波的真實波速。

2.4 數值各向異性

定義數值損耗各向異性Aα和數值相位各向異性Aβ分別為

式中：下標d表示數值波沿網格對角線方向傳播；下標a表示數值波沿坐標軸方向傳播。

在二維和三維情況下，PITD的數值損耗各向異性Aα和數值相位各向異性Aβ分別如圖4和圖5所示。電磁波頻率f＝300MHz，Courant常數S＝1 024，預先選定的整數N＝20，空間采樣密度Nλ＝20，電導率σ∈（1E－4，30）S／m，2D表示二維情況，3D表示三維情況。

在圖4中，三維情況下的Aα大于二維情況下的Aα，即PITD的數值損耗各向異性在三維情況下的值要大于其在二維情況下的值。在圖5中，三維情況下的Aβ的絕對值大于二維情況下Aβ的絕對值。當σ＝27.97mS／m時，Aβ＝0.PITD幾乎沒有數值波速的各向異性。在有耗介質中，FDTD的數值色散研究中也發現了類似的現象［15］。

綜合上述分析，有耗介質中PITD的數值損耗和數值色散都基本上不受時間步長的影響。隨著空間步長的減小，PITD的數值損耗和數值波速將趨近于電磁波的真實損耗和真實波速。當電導率較大時（σ≥ωε），PITD的數值色散特性要優于FDTD的數值色散特性。隨著維數的增加，PITD的數值損耗和數值色散的各向異性將增大。

3.數值算例

計算當繞組中電流突變時，變壓器鐵心疊片中磁場的瞬態變化。近似認為所有疊片中的磁場相同，只分析其中的一片。假設鐵心疊片的長度和寬度遠大于厚度，厚度方向為x方向，即磁場H是時間t和空間坐標x的函數。

鐵心疊片厚度a＝1mm，電導率σ＝1E7S／m，磁導率μ＝1.2E－3×πH／m.邊界條件

和初始條件

計算t＝5E－3 s時，在x＝0.25mm和x＝0.5mm處磁場強度H 的值。空間步長Δx＝a／1 024.在PITD中Courant常數S＝5E7，預先選定的整數N＝20.在FDTD中Courant常數S＝1.計算結果和CPU時間如表1所示。

表1 PITD和FDTD的計算結果和CPU時間

上述的計算結果表明，PITD的計算誤差要小于FDTD的計算誤差。相對于FDTD而言，PITD有更快的計算速度和更高的計算精度。由于采用了矩陣計算形式，PITD的計算內存需求較大。應用子域技術可以將PITD的內存需求降低至與FDTD相同的等級［16］。

5.結 論

有耗介質中PITD的數值損耗和數值色散特性的分析結果表明：PITD的數值損耗大于電磁波的真實損耗，其數值波速可以大于電磁波的真實波速。由于PITD色散方程中的時間量（τ）被限制在很小的范圍內，PITD的數值損耗和數值色散都基本上不受時間步長的影響。隨著空間步長趨近于0，有耗介質中PITD的數值色散方程將趨近于有耗介質中電磁波傳播的色散方程。因此，隨著空間步長的減小，PITD的數值損耗誤差和數值色散誤差都將逐步減小。當電導率較小時（σ?ωε），PITD的數值損耗和數值色散誤差大于FDTD的數值損耗和數值色散誤差。當電導率較大時，（σ≥ωε）PITD的數值損耗和FDTD的數值損耗趨于一致。由于在電導率較大時，PITD的數值損耗角正切與有耗介質的真實損耗角正切幾乎相同，使得PITD的數值波速比FDTD的數值波速更加接近于電磁波的真實波速。數值算例表明：對良導體而言，PITD比FDTD擁有更高的計算精度和更快的計算速度。

［1］葛德彪，閆玉波.電磁波時域有限差分方法［M］.2版.西安：西安電子科技大學出版社，2005.GE Debiao，YAN Yubo.Finite－Difference Time－Domain Method for Electromagnetic Waves［M］.2nd ed.Xi’an：Xidian University Press，2005.（in Chinese）

［2］TAFLOVE A，HAGNESS S C.Computational Electrodynamics：The Finite Difference Time Domain Method ［M］.3rd ed.Norwood：Artech House，2005.

［3］王秉中.計算電磁學［M］.北京：科學出版社，2002.

［4］李太全，田 茂，吳慶麟，等.圓環天線的時域有限差分分析［J］.電波科學學報，2002，17（6）：577－580.LI Taiquan，TIAN Mao，WU Qinglin，et al.Analysis of a round strip antenna using the finite difference time domain method［J］.Chinese Journal of Radio Science，2002，17（6）：577－580.（in Chinese）

［5］艾寶強，褚慶昕，雷振亞.基于FDTD算法的有源集成天線設計［J］.電波科學學報，2004，19（6）：739－741.AI Baoqiang，CHU Qingxin，LEI Zhenya.Design of active integrated antenna based on FDTD method［J］.Chinese Journal of Radio Science，2004，19（6）：739－741.（in Chinese）

［6］胡慧琳，譚云華，朱柏承，等.不同環境影響下共形蝶形天線性能的理論分析［J］.電波科學學報，2010，25（6）：1163－1168.HU Huilin，TAN Yunhua，ZHU Bocheng，et al.Conformal bow－tie antennas in different environments［J］.Chinese Journal of Radio Science，2010，25（6）：1163－1168.（in Chinese）

［7］SULLIVAN D M.Electromagnetic Simulation Using the FDTD Method［M］.New York：IEEE Press，2000.

［8］NAMIKI T，ITO K.3－D ADI－FDTD method unconditionally stable time－domain algorithm for solving full vector Maxwell’s equations［J］.IEEE Transactions on Microwave Theory and Techniques，2000，48（10）：1743－1748.

［9］FU Weiming，TAN E L.Stability and dispersion analysis for ADI－FDTD method in lossy media［J］.IEEE Transactions on Antennas and Propagation，2007，55（4）：1095－1102.

［10］MA Xikui，ZHAO Xintai，ZHAO Yanzhen.A 3－D precise integration time－domain method without the restraints of the Courant－Friedrich－Levy stability condition for the numerical solution of Maxwell’s equations［J］.IEEE Transactions on Microwave Theory and Techniques，2006，54（7）：3026－3037.

［11］ZHONG Wanxie，WILLIAMS F W.A precise timestep integration method［J］.Proc Inst Mech Eng，1994，208（6）：427－430.

［12］JIANG Lele，CHEN Zhizhang，MAO Junfa.On the numerical stability of the precise integration time－domain （PITD）method［J］.IEEE Microwave and Wireless Components Letters，2007，17（7）：471－473.

［13］CHEN Zhizhang，JIANG Lele，MAO Junfa.Numerical dispersion characteristics of the three－dimensional precise integration time－domain method［C］／／ Microwave Symposium IEEE／MTT－S International.Honolulu，2007：1971－1974.

［14］KRUMPHOLZ M，RUSSER P.On the dispersion in TLM and FDTD［J］.IEEE Transactions on Microwave Theory and Techniques，1994，42（7）：1275－1279.

［15］SUN Guilin，TRUEMAN C W.Numerical dispersion and numerical loss in explicit finite－difference timedomain methods in lossy media［J］.IEEE Transactions on Antennas and Propagation，2005，53（11）：3684－3690.

［16］白仲明，趙彥珍，馬西奎.子域精細積分方法在求解Maxwell方程組中的應用分析［J］.電工技術學報，2010，25（4）：1－9.BAI Zhongming，ZHAO Yanzhen，MA Xikui.Analysis and application of sub－domain precise integration method for solving Maxwell’s equations［J］.Transactions of China Electrotechnical Society，2010，25（4）：1－9.（in Chinese）