共形HIE-FDTD 方法的稳定性分析及改进

戴婧怡，杨 阳，李强兵，鲍卓如

(南京航空航天大学电子信息工程学院，江苏南京210016)

0 引言

1966年由K S Yee［1］提出的传统时域有限差分方法，已广泛应用于各种电磁仿真问题。为了提高计算速度、改善仿真效果和节省计算机资源，近年来，国内外学者提出了各种算法优化策略。然而FDTD方法本身存在着2个缺陷:①阶梯法模拟不能与网格共形的目标时，不可避免要进行阶梯近似，从而产生误差并可能会引起虚拟表面波;②柯朗－弗里德里希斯列维 (Courant－Friedrichs－Lewy，CFL)［2－3］稳定条件限制了时间步长的选取，使得计算效率受到限制。

对于阶梯近似的问题，Dey和Mittra［4］提出的金属共形时域有限差分方法，可以很好地改善由于阶梯近似产生的误差。但是，该方法并非没有弊端，在运用该方法时，需要满足2个稳定性条件，而且即便在满足条件的前提下，也不能做到时间步长的选取完全达到CFL条件。可以说该方法是以增加计算时间来换取较高的精度。对于CFL条件限制的改进，可以采用无条件稳定或弱条件稳定的方法，其中以交替方向隐式时域有限差分方法(alternating－direction implicit FDTD，ADI－FDTD)［5］和Crank－Nicolson时域有限差分方法(Crank－Nicolson FDTD，CNFDTD)［6－7］为主。1999年T Namiki提出的ADI－FDTD方法，将传统的时域有限差分方法中电场显式的迭代改成隐式迭代，这样时间步长的选取就不再受到CFL的限制，但是该方法的计算精度会随着时间步长的增加而降低。CN－FDTD方法在时间步长增大时离散误差较小，精度较高，但该方法需要求解一个大型稀疏矩阵方程组，在求解微细结构时效率很低。为了提高计算精度和效率，基于隐格式思想的 HIE－FDTD 方法［8－9］被提出来。该方法在电场迭代时，任意2个方向均采用隐式，另外一个方向采用显式，这样时间步长的取值可以与显式方向无关，即弱条件稳定。HIE－FDTD方法比ADI－FDTD方法有更高的精度，特别适用于在某个方向有微细剖分的结构。本文提出的金属共形HIE－FDTD方法采用x，y方向的隐式差分，着重分析在共形条件下减小HIE－FDTD时间步长的因子［10］，并给出修正方法。

1 Dey-Mittra的金属共形方法应用在HIE-FDTD

HIE－FDTD方法在1999年被提出，是基于无条件稳定思想的改进算法。算法采用电场分量在任意2个方向进行隐式差分，磁场分量差分格式不变的方法，可将CFL条件放宽，在第2节中进行推导［11］。

(1)—(2)式中:c为真空中光速;d为维数(对于二维问题，d=2;对于三维问题，d=3)。

本文分析的HIE－FDTD方法在x，y方向采用隐格式，z方向采用显式的格式。下面给出x方向和z方向E，H的离散方程为

同理，可以写出y方向电场与磁场的迭代公式。综上，采用z方向显格式x，y方向隐格式的HIE－FDTD方法的迭代公式，可以写成矩阵(7)的形式

当用HIE－FDTD方法模拟不规则金属边缘时，需要模拟的结构边缘不能与网格线重合，含有不规则边缘的金属板用阶梯网格划分边缘的截面示意图如图1a所示，其中，深色部分为金属，浅色部分为介质，syz为共形网格中介质部分的面积，ly，lz为共形网格中介质部分的边长。

图1 含有不规则边缘的金属片用阶梯网格划分边缘的截面示意图Fig.1 Cross－sectional view of a PEC boundary and the distorted cells

用Dey－Mittra提出的金属共形算法只需要对磁场表达式做如下修正。若含不规则边缘的截面只在某一个平面，则只有一个方向的磁场需要修正，以x方向为例，磁场表达式由(5)修正为

Dey－Mittra的金属共形方法需要满足2个条件以保持算法的基本稳定性［12］。

1)共形网格面积s(i，j，k)应大于规则网格面积δ2的5%。

2)共形网格中的最长相对边长l'(i，j，k)与该网格相对面积 s'(i，j，k)的比率小于12。

当满足上述条件，Δt也需要适当减小，以保证稳定。因此，将该共形方法应用在HIE－FDTD方法中，也会导致时间步长的减小，或发散的提前。下面就导致发散提前的不稳定因子进行分析，并给出改进方法，进行数值验证。

2 金属共形HIE-FDTD方法的稳定性分析

HIE－FDTD方法的时间步长的选取只与隐格式方向的空间步长有关，下面进行简要的推导。

参考HIE－FDTD方法的时域基本推进公式的矩阵形式(7)，考虑平面波本征模的解以及有限差分近似，即

(9)式中:f(i，j，k)=e0xyz;ku为波矢量;ζ为增长因子;u=x，y，z;V=E，H 带入矩阵(7)式得到新的矩阵形式为

(11)式中，ru=(cΔt)2W2u。

解得

为了保证数值结果的稳定性，要求增长因子|ζi|≤1，即

经过推导可以得出:HIE－FDTD方法的数值稳定性条件为

将变形后的修正磁场表达式(15)带入矩阵(7)得到

为使该矩阵方程有非零解，其系数行列式必须为零，因此，我们得到与增长因子ζ有关的多项式为

K是由同一个网格中的介质与金属的分界产生的，K的大小取决于小边长与小面积的比值。当网格为规则网格时，K≡1，(16)式即还原为常规HIEFDTD方程，即时间步长的选取只需满足(14)式。

Dey－Mittra的方法中，在图1b的情况下，有

在图1c的情况下，有即，在Dey－Mittra方法中，需要分情况讨论且不能保证稳定性。要使共形方法尽量不破坏原迭代方程的稳定性，则需要因子K趋近于1。

根据K的计算公式，要使得K≡1，可令

则

容易求得:syz(Dey－Mittra)≤sequ。更容易满足稳定条件1)。根据(20)式，可以看出本文提出的方法较Dey－Mittra的方法可以保持HIE－FDTD固有的稳定条件。

3 数值证明

一个功分器结构［8］如图2所示。器件尺寸分别为 W50=4.5 mm，L50=10 mm，d=4.3 mm，hl=4.73 mm，Wl=28.4 mm，ll=11.92 mm，介质板厚度为1.575 mm，相对介电常数 εr=2.2。取dx=0.787 5 mm，dy=1.5 mm，dz=0.1 mm。图3为端口2的时域波形，图3上的曲线分别是应用Dey－Mittra 的方法分别取 dt=2.311 3 ps=7ΔtHIE－FDTD，dt=1.650 9 ps=5ΔtHIE－FDTD，dt=1.320 7 ps=4ΔtHIE－FDTD，以及应用本文提出的方法，取 dt=2.311 3 ps=7ΔtHIE－FDTD的端口时域波形，从图3中可以得出结论:①减小时间步长可有效地推迟发散;②本文提出的方法，较Dey－Mittra的方法，可以大大推迟发散，在时间步长取相同的情况下，可将发散出现的时间步推迟3.18倍，(即，使时间步长取 Dey－Mittra方法的1.75倍，发散可推迟1.25倍，时间步长取Dey－Mittra方法的1.4倍，发散可推迟1.4倍)。图4给出了本文提出的算法在取 dt=2.311 3 ps=7ΔtHIE－FDTD以及 Dey－Mittra 算法取 dt=1.320 7 ps=4ΔtHIE－FDTD时端口2的S21和仿真结果相比较的对比图。可以看出，本文提出的方法具有较高的精度。表1给出了本文提出的算法取dt=2.311 3 ps=7ΔtHIE－FDTD以及 Dey－Mittra 算法取 dt=1.320 7 ps=4ΔtHIE－FDTD时，CPU time 的比较，从表1 可以看出，时间节约878 s(约42%)。结果证明:本文提出算法，能有效推迟发散，并具有较高的精度和效率。

表1 Dey－Mittra 方法取 dt=1.320 7 ps=4ΔtHIE－FDTD时和本文提出方法取 dt=2.311 3 ps=7ΔtHIE－FDTD时的 CPU 用时比较Tab.1 CPU time of Dey－Mittra when dt=1.320 7 ps=4ΔtHIE－FDTD and proposed method when dt=2.311 3 ps=7ΔtHIE－FDTD

图2 功分器结构Fig.2 Geometry of the conventional TPD

图3 端口2的时域波形图Fig.3 Wave shape of time domain electric field

4 结论

本文将Dey－Mittra的金属共形方法应用在HIEFDTD，并分析了该金属共形HIE－FDTD方法不能取到CFL条件的原因，提出了一个可以增加算法稳定性的面积计算方法。数值实验结果证明:提出的方法易于实现，当2种方法取相同时间步长，本文提出方法较Dey－Mittra的算法发散出现时刻推迟3.18倍，在相同精度情况下，计算时间可以节约42%，是一种高效的计算方法。

图4 Dey－Mittra共形算法、本文提出算法和文献结果的对比图Fig.4 Comparation of S21with the Dey－Mittra’s method，the proposed method and the reference result

［1］YEE Kane S.Numerical solution of initial boundary value problems involving Maxwell＇s equations in isotropic media［J］.IEEE Transaction on Antennas and Propagation，1966，14(3):302－307.

［2］TAFLOVE A，HAQNESS S C.Computational Electromagnetics:The finite－difference time－domain method ［M］.Boston，MA:Artech House，2000.

［3］童创明.计算电磁学快速方法［M］.西安:西北工业大学出版社，2010.TONG Chuangming.Fast method for computational electromagnetics［M］.Xi＇an:Northwest Polytechnic University Press，2010.

［4］DEY Supriyo ，MITTRA Raj.A locally conformal finite difference time－domain(FDTD)algorithm for modeling three－dimensional perfectly conducting objects［J］.IEEE Microwave and Guided Wave Letters，1997，7(9):273－275.

［5］NAMIKI Takefumi.A new FDTD algorithm based on alternating－direction implicit method ［J］.IEEE Transactions on Microwave theory and technology，1999，47(10):2003－2007.

［6］YANG Y，CHEN R S，WANG D X，et al.Unconditionally stable Crank－Nicolson FDTD method for simulation of three－dimensional microwave circuits ［J］.IET Microw.Antennas Propag，2007，1(4):937－942.

［7］YANG Y，FAN Z H，DING D Z，et al.Extend two－step preconditioning technique for the Crank－Nicolson finitedifference time－domain method to analyze the 3D microwave circuits［J］.Int J RF Microw C E，2009，19(4):460－469.

［8］van VIEN，CHAUDHURI Sujeet K.A hybrid implicit－ex－plicit FDTD scheme for nonlinear optical waveguide modeling［J］.IEEE Transactions on Microwave Theory and Techniques，1999，47(5):540－545.

［9］HUANG B K，WANG G，JIANG Y S，et al.A Hybrid Implicit－Explicit FDTD scheme with weakly conditional stability［J］.Microwave and Optical Technology Letters，2003，39(2):97－101.

［10］DAI Jian，CHEN Zhizhang，SU Donglin，et al.Stability Analysis and Improvement of the Conformal ADI－FDTD Methods［J］.IEEE Transactions on Antennas and Propagation，2011，59(6):2248－2258.

［11］兰婧，杨阳，戴婧怡.FDTD全波分析缺陷地结构微带线的高效算法实现［J］.重庆邮电大学学报:自然科学版，2012，24(4):426－430，466.LAN Jing，YANG yang，DAI Jingyi.Highly－efficient FDTD method of full－wave analysis on the microstrip defected ground structure［J］.Journal of Chongqing University of Posts and Telecommunications:Natural Science Edition，2012，24(4):426－430，466.

［12］葛德彪，闫玉波.时域有限差分方法［M］.2版.西安:西安电子科技大学出版社，2002:184－185.GE Debiao，YAN Yubo.Finite－Difference Time－Domain Method for Electromagnetic Waves［M］.2nd ed.Xi＇an:Xi dian Press，2006:184－185.