《计算电磁学》第六讲

1、0:52:07计算电磁学第六讲Dr. Ping DU (杜平)E-mail: School of Electronic Science and Applied Physics， Hefei University of Technology (HFUT)Nov. 3, 2011FDTD中的若干技术及应用0:52:07 OutlineI. 激励源的设置II. 集总元件的模拟III. 近场-远场变换IV.FDTD法的一些应用实例0:52:07I. 激励源的设置为减少引入源带来的内存占用和计算时间，要求激励源的实现要紧凑.1.1 强迫激励源 在FDTD网格中，通过直接对特定的电场或磁场分量强行赋予所需

2、的时间变化形式，可简便地建立强迫激励源。 比如，对一维TM网格，波沿x轴传播，电场只有 分量，为模拟频率为zE、在n=0时开始的连续正弦波源， 0f可在源处( )令 为 sizE00( )sin(2)nzsEiEf n t另一种常用的源是高斯脉冲 220000( )expexpnzsdecayttnnEiEETn (6-1)(6-2)0:52:07若要求高斯脉冲在t=0时刻近似为零，则应选取 03tTT由所需的脉冲频谱带宽决定 高斯脉冲的频谱 222( ) expG fT f也是高斯形的。超过某个频率后，其频谱的幅度很小了，该频率称为 maxf如，可选 max1=2fT如果要使直流分量为零，有

3、效中心频谱中心位于 0f,则可采用调制高斯脉冲作为源20000( )expsin 2()nzsdecaynnEiEfnntn(6-4)(6-3)(6-5)0:52:07上述强迫激励源会产生相应时间变化的数值波, 将从源点 si向两个方向对称传播 于是得到了形式紧凑的激励源。 抽样函数激励源sin(/ 2)( )/ 2tg tt其傅立叶变换为 2( )44G fUfUf频谱主要分布在主峰附近，远离它的频谱贡献可忽略不计，可只取主峰前后2T范围内的波形。考虑到FDTD执行时间从0开始，可将时域波形平移为 (6-7)(6-6)0:52:07sin()/ 2( )()/ 2tTg ttT当研究单模工作

4、波导等具有一定工作频带电磁结构时，该频谱限带信号最合适 设波导单模工作带宽为 2，中心频率为 cf，则可取 0sin()/ 2( )cos(2)()/ 2nzsctTEiEf ttT(6-8)(6-9)0:52:07表7-1 稳态源激励下场达到稳定所需的时间散射体结构类型稳定所需波源周期数凸型二维金属散射体，跨度小于1,TM波有耗三维结构，特别是由生物体组成的三维结构5凸型二维金属散射体，跨度为15凸型三维散射体，跨度为15520三维金属线或棒，跨度为1，接近谐振激励三维金属散射体，跨度为10，有角反射和开放腔2040深度重入式三维金属散射体，跨度1040三维任意尺寸金属散射体，具有中等到高Q

5、的带孔谐振腔，接近谐振激励1000:52:07 为克服强迫激励源引起虚假反射, 在激励脉冲几乎衰减为零、来自结构的反射波到达 网格点之前，将激励去掉。 网格点场值的刷新换用标准FDTD公式。如果 网格点为边界点，场值的刷新换用吸收边界条件。 sisi但这一办法，对(6-1)的正弦激励源不合适.因为激励在来自结构的反射波到达 网格点之后仍然存在，不能被去掉. si另外，对脉冲宽度较宽的高斯脉冲(包括调值和非调制)激励源，为使来自结构的反射波到达 网格点之前激励脉冲几乎衰减为零,网格点到结构的距离要充分大，这会增加内存占用、增加计算时间。 si 但由于形式简单、使用方便，强迫激励源还是被用来处理许

6、多工程问题，只要其使用条件满足. 如波导系统中同轴线探针的模拟，就可将强迫激励源放置于金属探针处。 0:52:07克服强迫激励源虚假反射的另一办法: 将激励源看成有源Maxwell方程中的电流密度J进行处理。根据Maxwell方程11stEHJ分量的FDTD差分格式为 zE11212,( )( )()( )snnnnzszszz ssittEiEiJiH记等效电场激励源为 112,( )( )nnz ssz sstEiJi 总场FDTD公式为 1211,( )( )()( )snnnnzszszz ssitEiEiEiH这种引入激励源的方式不会产生虚假反射(6-10)(6-11)(6-12)(

7、6-13)0:52:07注意：等效激励源 1,( )nz ssEi并不等于 szi处的电场值 1( )nzsEi，它仅仅是激励 1( )nzsEi的源.由于这种激励以FDTD公式为载体，因而不能放置在仿真区的边界。这里的附加激励源是在电场FDTD中引入，相当于电压源激励；如果在磁场FDTD中引入，相当于电流环激励。 强迫激励源常被用来模拟波导系统中的入射波。 根据所需模拟的传播模式的时空分布，取某一波导横截面为激励源平面，对位于该面上的电场切向分量赋予相应的时空变化。此时，要将(1)、(2)、(5)中的替换为相应模式的横截面的场分布函数。 如采用高斯脉冲激励，波导系统的传输方向沿z轴，模式的横

8、截面场分布为 tan( , )nEx y, 强迫激励源加在 处 szkz于是，强迫激励源为 0:52:072200tan0,tan0,tan( , ,)( , )exp( , )expnsdecayttnnEi j kEi jEi jTn这样方式不会激起不希望存在的模式。 若模式的横截面场分布未知，有两个办法： (1)采用近似的模式分布； (2)通过预处理，数值仿真获得所需模式分布。 可能激励不希望存在的数值模式。为此需要增加一长段辅助波导让高次模充分衰落，余下较纯的所需模式分布。 先用FDTD仿真一个与要研究的波导截面相同、足够长的均匀波导，使高次模充分衰落；然后在波导远端记录横截面切向电场

9、分量的分布。从而得到了强迫激励源。 需要注意：激励源的频谱尽量不包含低于截止频率的分量。 (6-14)0:52:071.2 总场/散射场体系 在FDTD中，它用的最多。以电磁散射问题为例进行说明，如图1. 图1 将FDTD网格划分为总场区和散射场区 0:52:08 将仿真区域分为1区和2区，前者为总场区(入射场+散射场)，在该区用总场FDTD公式；后者为散射场区，用散射场FDTD公式。区域2的外边界为截断边界，用吸收边界条件吸收外向波。连接边界上的场为总场。 只要入射场的时空变化准确给出，上面的激励源引入过程不含任何近似，也不引入任何虚假模式，且形式紧凑，内存占用少。 例1 考虑金属机箱，左面

10、开孔，位于TEM平板传输线中，TEM平面电磁波照射到该机箱，一部分波通过孔进入机箱，其余被反射，如图2所示。 图2 机箱屏蔽问题 0:52:08设入射波只有 、 分量 ,in xE,0/in yin xHE激励连接面放在距离机箱含孔面为L的地方，并以此平面将场区分为散射场区(不含机箱)和总场区(含机箱)。 为保证电磁场的连续性，在此平面需引用激励连接条件。 设激励连接面位于 处、 和 位于激励连接面上且取为总场。 zk z xEyE激励连接条件为：在连接平面上 112,1122,1122,1111, , , ,22221111,22221111, , ,2222nnnt xt xin ynnt

11、 zt znns yt ytEij kEij kHij kzHijkHijktyHij kHij kz(6-15)0:52:08112,1122,1122,1111,22221111,22221111,2222nnnt yt yin xnnt xs xnnt zt ztEi jkEi jkHi jkzHi jkHi jktzHijkHijkx(6-16)0:52:08在连接面左邻半个步长的网格点上，112,11111,2222211,12211, ,1,22nnns xs xin ynnt ys ynns zs ztHi jkHi jkEi jkzEi jkEi jktzEi j kEi jk

12、y(6-17)0:52:081122,11111, , , ,222221111, ,22211, ,1, ,22nnns ys yin xnns zs znns xt xtHij kHij kEij kzEijkEi j ktzEij kEij kz 其中，t、s、in分别表示总场、散射场、入射场。 本例中， ， ,0in xH,0in yE(6-18)0:52:08在总场/散射场体系中，激励源的引入很紧凑，只要知道 0zkz平面上的切向电场和012zkz平面的切向入射磁场。 对一般的导波系统，设入射波为高斯脉冲，模式横截面电场分布为 0,tan( , )x yE、横截面磁场分布为 0,ta

13、n( , )x yH，激励源可表示为 0200,tan00,tan200,tan, ,( , )exp( , )expninz zdecayzzttvi j ki jTnni jnEEE(6-19)0:52:08020102,tan00,tan200,tan12, ,( , )exp2122( , )expninz zdecayzztttvi j ki jTznnv ti jnHHH其中，v是波模沿z轴的传播速度。 对无色散的波模， 1v对无色散(或弱色散)的TEM(或准TEM)传输线中的主模，若其等效介电常数 eff已知，则 . 1effv对一般波模，v可以事先通过一个FDTD数值模拟过程计

14、算获得。 (6-20)0:52:08总场/散射场体系具有如下特点： 可实现任意波激励可直接使用吸收边界条件易于计算远场响应更宽的计算动态范围不适合于分析宽带色散导波系统 对色散严重的导波系统， 是频率的函数，式(6-19)中的v取值不确定. ( )v0:52:082. 集中参数电路元件的模拟 tDLJ在Maxwell方程中的位移电流之后再引入一项集总电流密度可以计及集总电路元件的影响。方程扩展为 ,LtDHJ, ,zi j kELJLI假设元件位于自由空间，沿z轴放置，与电场元件的截面积按一个网格计，则元件的局部电流密度和元件总电流满足 处于同一位置。LLIJx y 2.1 扩展FDTD方程,

15、 xy其中分别为x、y方向的步长， LI是元件两端电位差 , ,Vzi j kEz的函数 规定的 正向为+z方向。 LI(6-21)(6-22)0:52:08于是得到扩展FDTD方程 11212, , ,00, ,()nnnnzzzLi j ki j ki j kttEEIx y H对于沿x轴和y轴取向的集总参数元件可类似处理。 2.2 几种集总参数电路元件 (1) 电阻 用FDTD分析电路时，可能会遇到含电阻的情况。可以采用实体电阻；也可以采用前述的集总参数元件模型。 对于电阻，其电流-电压特性为 11, , ,2, ,()2nnzzni j ki j kzi j kEEzIR(6-23)(

16、6-24)0:52:08将其代入式(6-23)，整理得 1, ,nzi j kE的时间步进关系 11002, , , ,0012()1122nnnzzzi j ki j ki j kt ztRx yEEt zt zRx yRx y H 为验证该集总参数电阻在FDTD计算中的性能，有人曾模拟了50微带线端接匹配电阻的例子，激励源是一有效带宽为20GHz的高斯脉冲。数值实验表明，直到1GHz范围内，这两种电阻模型产生的反射系数均小于1%. (6-25)(6-26)12, ,nzi j kLIJx y 0:52:08(2)阻性电压源 利用扩展FDTD方程模拟集总参数元件的能力，可以很方便地模拟一个无

17、反射的阻性电压源。 阻性电压源的电流-电压特性为 1112, , ,2, ,()2nnnzzni j ki j kszi j kssEEzVIRR12, ,nzi j kLIJx y 其中，12nsV是电压源 ，sR是源内阻。 将其代入式(6-23)，整理得 1, ,nzi j kE的时间步进关系 (6-27)(6-28)0:52:0811122, , , ,01 AB()1A(1)1B/ 2nnnnzzzsi j ki j ki j ktEEVAzH其中， 0A2t zRx y ，0BstRx y (3)电容 电容的电流-电压特性为 11, , ,2, ,()nnzzni j ki j kz

18、i j kEEzICt12, ,nzi j kLIJx y 其中，C是电容值。 (6-29)(6-30)(6-31)0:52:08将其代入式(6-23)，整理得 1, ,nzi j kE的时间步进关系 1102, , , ,0()1nnnzzzi j ki j ki j ktEEC zx y H(4)电感电感的电流公式为 0tVIdtL其中，L为电感值。 电感的电流-电压特性为 (6-32)(6-33)0:52:0812, , ,1nnmzzi j ki j kmz tIEL 12, ,nzi j kLIJx y 将其代入式(6-23)，整理得 1, ,nzi j kE的时间步进关系 2112

19、, , , , ,100()()nnnnmzzzzi j ki j ki j ki j kmtztEEELx y H(5) 二极管集总参数二极管的电流-电压特性为 0exp1ddqVIIkT(6-34)(6-35)(6-36)(6-37)0:52:08其中，q是电子所带的电量，Vd是二极管两端电压，k是玻尔兹曼常数，T是绝对温度。 将其代入式(6-23)，整理得 1, ,nzi j kE的时间步进关系 112, , , ,010, , ,0()exp()12nnnzzzi j ki j ki j knnzzi j ki j ktEEtq zIEEx ykT H 利用牛顿法求解（6-38）可得电

20、场刷新值。数值实验表明，该模型在直到15V的二极管电压范围内都是数值稳定的。 (6-38)0:52:08(6)双极型晶体管 对NPN双极型晶体管，也能用FDTD算法实现大信号分析，包括对数字开关过程的分析。图3所示为一接在微带线上的双极型晶体管，发射极接地，基极接在微带线上. 图3 端接在微带线上的共发NPN双极型晶体管的三维FDTD模型 0:52:08根据Ebers-Moll晶体管模型，其电流-电压满足 0exp1BEFqVIIkT0exp1BCRqVIIkTERFIIICRFFIII 如图3所示，假设晶体管位于自由空间，沿z轴取向，微带与接地板的间距是一个空间步长。 BEVzEBE基极-发

21、射极电压可用FDTD法算得的电场来表示 (6-39)(6-40)(6-41)(6-42)0:52:08112()2nnnBEzzEBEBzVEE BCVzBCE相似地，基极-集电极电压可用FDTD法算得的电场来表示 112()2nnnBCzzBCBCzVEE将式(6-43)、(6-44)代入到式(6-39)-(6-42)，可得 1112010exp12exp12nnnERzzBCBCnnzzEBEBq zIIEEkTq zIEEkT(6-43)(6-44)(6-45)0:52:081112010exp12exp12nnnCzzBCBCnnFzzEBEBq zIIEEkTq zIEEkT将这两式

22、代入式（6-23），可得晶体管处电场分量的FDTD刷新公式 1112200()nnnnzzzEEBEBEBttEEIx y H1112200()nnnnzzzCBCBCBCttEEIx y H用牛顿-拉夫森法可求得上述耦合超越方程组。数值实验表明，用上述方法算得的与用SPICE模型求得的解吻合很好。 牛顿-拉夫森方法(Newton-Raphson method)： 用牛顿-拉夫森方法求解 ( )0f x 的根。迭代公式为 (6-46)(6-47)(6-48)0:52:081(),0,1,2,()nnnnf xxxnfx使用时，要先设定一初值 。 3. 近场-远场变换 (Near-Field t

23、o Far-Field Transformation) 0 x 采用总场/散射长体系，散射场区的外边界可以不随散射或辐射体而变化，因而可以用统一不变的计算过程由外边界上的近区场计算出远区散射或辐射场。先考虑时谐场(time harmonic field)情形。3.1 二维TM波 场分量有3个：zxyEHH、设FDTD仿真区域(散射场区)的外边界为C，nc为C的单位外法向矢量， C上的切向电场和切向磁场分别为 tantan,CCEH，“”表示频域量。 0:52:08根据等效原理，C上的切向等效电流和等效磁流分别为 tanCCeqCJnHtanCCeqC MnE频域远区场为 ( /4)0( )(

24、)( )8limjkjCCjkzzeqzeqkCeeEkedCkaaa JM其中， zaa、分别为沿z轴和径向的单位矢量。 如果定义复数方向图函数为 ( /4)0( )( )( )8jCCjkzeqzeqCeFkedCkaaa JM(6-49)(6-50)(6-51)(6-52)0:52:08则双站雷达散射截面(RCS)为 22,( )( )2z incFRCSE为计算远场，需先知道C上的 tantan,CCEH。若切向电场位于网格外边界C上，等效磁流可由C上的 tanCE算出； 等效电流需要 tanCH，其可以从C后退半个步长的网格点上取值。 为获取频域数据，可以利用离散傅立叶变换。 ()G

25、 n t设表示网格点上切向场的时域采样。其频域值为 102()()exp, (0,1,2,)NfnjknG k ftG n tkNN 其中N为总的步进数，f=1/(Nt)。 (6-53)(6-54)0:52:083.3 三维情形 设FDTD仿真域的外边界是S，其单位外法向矢量为ns， S上的切向电场和切向磁场分别为 tantan,SSEH由等效原理，S上的切向等效电流和等效磁流分别为 tanSSeqSxxyyzzJJJaaaJnHtanSSeqSxxyyzzMMM aaaMnE频域远区场可由下式计算 ( , , )0rE r 0( , , )()4jkrjkeE rLNr 0( , , )()

26、4jkrjkeErLNr (6-55)(6-56)(6-57)(6-58)(6-59)0:52:08( , , )0rHr 0( , , )()4jkrLjkeHrNr 0( , , )()4jkrLjkeHrNr 其中， 000/为自由空间波阻抗， NNLL、 、按下列式子计算 cos( , , )coscos( , , )cos sin( , , )sinxyjkrSzJx y zJx y zNedSJx y zcos( , , )sin( , , )cosjkrxySNJx y zJx y zedS(6-60)(6-61)(6-62)(6-63)(6-64)0:52:08cos( , ,

27、 )coscos( , , )cos sin( , , )sinxyjkrSzMx y zMx y zLedSMx y zcos( , , )sin( , , )cosjkrxySLMx y zMx y zedS其中， 是远区场点位置矢量r和等效面元S上源点位置矢量 r之间的夹角. 假设S是一个边长为 的矩形盒子，中心在原点。则上述闭合面上的积分可氛围3组面积分： 000222xyz、 、(1) 位于 0 xx 的两个表面 SeqJ和SeqM的非零分量： yzyzJJMM、 、指数相位项： 。(6-65)(6-66)0:52:0800cos()(sincossinsincos )sincoss

28、insincosrxyzxyzrxyzxyz r aaaaaaa积分区间: 0000,;yyyzzz dSdydz。(2) 位于0 yy 的两个表面 SeqJ和SeqM的非零分量： yzxzJJMM、 、指数相位项： 00cos( )(sincossinsincos )sincossinsincosrxyzxyzrxyzxyzr aaaaaaa积分区间: 0000,;xxxzzz dSdxdz。(6-67)(6-68)0:52:08位于的两个表面 SeqJ和SeqM的非零分量： (3)0 zz yxyxJJMM、 、.指数相位项： 00cos( )(sincossinsincos )sinco

29、ssinsincosrxyzxyzrxyzxyzr aaaaaaa积分区间: 0000,;yyyxxx dSdxdy.散射场的时间平均Poynting矢量为 *2220022011Re()Re()22(| )32scatPE HE HkLNLNr (6-69)(6-70)0:52:08则双站RCS为 2222000( , )lim 4(| )8sactrincincPRCSrPkLNLNP 其中， 为入射波功率密度。对于辐射问题(如天线/天线阵)，式(6-71)仍然可以用来计算远区场方向图。因为此时没有入射波， 可以取incPincP天线的总入射功率。(6-71)0:52:083. FDTD法

30、应用例子 3.1 均匀三线互连系统3.2 波导元件的高效分析3.3 传输线问题的降维处理3.4 周期性结构分析 0:52:083.1 均匀三线互连系统 互连系统设计是目前超高速集成电路研制中的一项关键技术。图4所示为一均匀三线互连系统。当超宽频使用的时候，信号的有效频谱从DC一直到微波段。加上超微细结构，他的特性与低频电力传输线、高频微波传输线有同有异。 图4 均匀三导体互连系统 0:52:09 采用FDTD法分析多导体互连线的传输特性，将导体条带及接地板作为导电率为 的媒质处理，能计及趋肤效应、邻近效应以及色散等因素的综合作用。 本问题的FDTD求解过程为： 1)在初始(t=0)时令所有的场量都为零 . 2)将一高斯脉冲激励 加在端口1。 220exp() /yEttT先从FDTD方程算出 时刻的磁场强度 (1/ 2)tnt12nH再从FDTD方程算出 时刻的磁场强度 (1)tnt1nE令接地板下表面切向电场为零，并在截断边界上用吸收边界条件；0:52:09 记录