时域有限差分法姚伟介绍

1、伊犁师范学院硕士研究生期末考核科目：电磁波有限时域差分方法姓名：姚住学号：1076411203009学院：由子与信息工程学院专业：无线电物理时域有限差分法1选题背景在多种可用的数值方法中，时域有限差分法 （FDTD）是一种新近发展起来的可选方法。 1966年，K.S.Yee首次提出电磁场数值计算的新方法一时域有限差分法（FiNte Difference- Time Domain ,简称FDTD）经历了二十年的发展 FDTDt才逐渐走向成熟。 上世纪80年代后期以来FDTDt进入了一个新的发展阶段，即由成熟转为被广泛接受和应 用的阶段。FDTDt是解决复杂问题的有效方法之一，是一种直接基于时域电

2、磁场微分方程的数值算法，它直接在时域将 Maxwell旋度方程用二阶精度的中心差分近似，从而将时域 微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计 算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机 内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量 学、瞬态电磁场研究等多个领域1。2原理分析2.1 FDTD 的 Yee元胞E,H场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理.D 壬=S.:t ft'

3、、E =汨 ,:H=r1.:t 2 t图1 Yee模型如图1所示，Yee单元有以下特点2:1） E与H分量在空间交叉放置，相互垂直；每一坐标平面上的E分量四周由H分量环绕，H分量的四周由E分量环绕；场分量均与坐标轴方向一致。2）每一个Yee元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用 棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。3）每一场分量自身相距一个空间步长，E和H相距半个空间步长4）每一场分量自身相距一个时间步长，E和H相距半个时间步长，电场取n时刻的值， 磁场取n+0.5时刻的值；即：电场n时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值 由

4、n-0.5时刻的值得到；电场n时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的旋度对应(n+0.5)+0.5时刻的电场值，逐步外推。5) 3个空间方向上的时间步长相等，以保证均匀介质中场量的空间变量与时间变量完全对称。应用这种离散方式，将含时 间变量的Maxwell方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。 由电磁问题的初值和边界条件，就可以逐步推进地求解以后各时刻空间电磁场分布。2.2 Maxwell方程FDTD勺差分格式VxH ="+J麦克斯韦第一、二方程 t(1)' E = - - J mT式中，J时电流密度，反映电损耗，Jm是磁流密度，单

5、位V/m反映磁损耗。主要 与上式对应。各向同性介质中的本构关系：(2)D = E B = H J = EJm = m H其中力是磁阻率，计算磁损耗的以E，H为变量，在直角坐标中，展开麦克斯韦第一、二方程，分别为-：Ey什-：Ey-Ez:t:t-HxEyEz(3)二 y-Ex:zft- mHx-Ez-：Hy：z'EyEx-yft=_ J以-:t- mH y-'mHz令f(x,y,z,t淤表E, H在直角坐标中的任何一个分量，离散符号取为(5)f x,y,z,ti；=f i x, j y,k z,n ：t )= fn i,j,kf(X, y,z,t 次于时间和空间的一阶偏导数取中心

6、差分近似为 定;ifn(iT，j,k)_fn。：，"x=i4&X wf n(,j +,k)-f n(i,j -4，kjl(6),矽心 wf n(i , j ,k +； )- f n(i, j ,k 一： N友 z4& Az尸f1 n+n-1»- |f 2(i，j*f 2(i,j,k)、ct t 朝& |可以看出，每一节点上沿某一方向场分量的一阶偏微分可以用在该方向上相邻两点的 一阶中心差商来描述，将式（1）用一阶中心差商方程取代，整理后便得到一阶差分方程， 它具有二阶精度3 oYee元胞如图1所示，规定为1）剖分节点与场分量所在棱边中点不同，场分量的

7、位置，即 E,H节点是Yee元胞节 点的相对位置，不需要单独编码；2）当空间存在媒质分界面时，场量自动满足场的连续性条件，E1t = E2t，H1t = H2t电磁分量的取样方式不仅符合法拉第电磁感应定律和安培环路定律的自然结构，也符 合麦克斯韦方程的差分计算。其次，时间步长可以取为电磁波传播一个空间步长所需时间 的一半，因此E与H在时间顺序上交替抽样，时间间隔相差半个时间步长。2.3 一维问题均匀平面波（TEM波）是一维问题，设电磁波沿 z轴方向传播，则Ez=°，Hz = ° ,cc=0，一 =0场量和介质参数均与x, y无关，即敌到,麦克斯韦方程为.:E:£；

8、：t汩Uft-+ ymH y；H y _;z：EX 二 一z-：Hx:z 王y:z,:E£ft.:H p：(8)-:t旋转坐标轴后可以只保留一组公式4,设保留（7）Yee元胞如图2所示差分格式为ExHy图2 一维Yee元胞c1E；Rk)=CA(mE(k)CB(m)1 Hz1+1-h y2 k -1(9)11nn _Hy 2 k 1 : CP m Hy 2 k 2 -CQ m（i0）如果介质无损耗，则=0，m = 02.4二维问题三维通常是散射问题，二维是 TE、TM波问题，一维是TEM发问题。在二维场中，所有物理量与 Z坐标无关，既已/金=0。于是在TE和TM波的表达式分别为TE波(

9、Ez=0).:Hz：:Exz 二；一xEx.y .t.,Hz-Ey Eyx .t沱 y FEx -Hz改 ：:y：:t一 mHz(11)TM波(Hz =0)：Ez ：Hx一 mHx,t:H Vy - mHyFt)；：Ez Ezft:y-：Ez：x'Hy.x图3分别给出了 TM波和TE波的Yee元胞图图3 TM波的Yee元胞=3(12)：:H；:yEy(i+1/2, j+1/2)Hz (i, j+1/2)Ex图4 TE波的Yee元胞对于TE波，只要令巳二°,在Az上，Hx,Hy不随z变化，m中去掉k即可得至IJ:111、n a 1n 1nnt -inn七 1ExT +2J )=

10、CA(mEx(i+2，j )+CB(m) Hz 2(i+2,j + 2 )-H z 2(i+2,j-2)八J式中：m = i +2, j(13)/njn i、!-n 书,1 ca/ u-nf. .上 1 o k 1- - 0 2 f , 1 , 1 0 2 (1 , 1 Ey q, j + 2 )=CA(m Ey i, j + 2 )-CB(m )- H z 2Q + 2 , j + 2 ) H z 2 U _ 2, j + 2 )八)式中：m =i, j +：(14)11Hz 2(i +2, j + 2 )=CP(mHz 7 +2，j + 2)E；(i +1,j +； )-E；(i,j +；

11、 ) EX(i+： , j +1)-EX(i+1 , j ) -CQ(m 1AxAy式中，m=i +2 ,j +2(15)对TM波，只要令Hz =0 ,在Az上，Ex,Ey不随z变化，m中去掉k,即可得到:n Jn 1AhX 2 i,j ； =CPmH；W i,j .； -CQm - E； i,j 1 -E； i,j 1y式中,m=i, j +：(16)n/n 1H；消 +2，j )=CP(mH>(i +2，j)+CQ(m) 1 13+1/)E；(i,j )1x(17)式中，m =i jEL i,j =CAm E； i,j11111h；2G+.)h；2G：,j) H：2(i,j+1AHX

12、2(i,j(18)式中：m =i, j为了编写统一的TE和TM波二维FDTLg序，可将描述TE波差分公式(13)(15)中相 应的标号整体移动1/2 ,即坐标(x,y )分别沿x和y轴方向移动半个网格，并将离散时间 也移动半个时间步长，式(13)(15)可以重新写为n 1n1Ex 2 i,j J =CAmEx 2 i,j 弓 CB m H； i,j 1 -H； i,j 1(19)11n1nEy 2 i 2-,j 二CA m Ey 2 ijCBmini 1'j -Hl" 1(20)Hn1 i,j =CPmHn i,j-CQ(m)1n 5Ey 2 in 51n .j -Ey 2

13、i -7,jEx 2 i,jx1n2 1_ Ex i,j - 2y式中：m =i +；, j(21)式中，m=ij可以看出，TE波的FDT必式(19)(21)与TM波的FDT必式(16)(18)形式相同, 给编程带来极大方便。注意TE波和TM波之间的对偶关系5 ,即这样就可以编写统一的计算程序了2.5三维问题(直角坐标系)2.5.1 电场时间推进差分格式(i +：, j,k )、(i, j +k X (i,j,k 十;位置上的 Ex,Ey,Ez节点(i,j,k的3个电场分量分别用 表示，以式(3)中第一个公式为例：必也:.y；z在"*"+："问步，对节点(i,j,

14、k )的离散公式为:ei 1 ,j,k -J i ； ,j,k -E'i ； ,j,k；,j,k - tn 工in 乂n i J jkEnijkH 2 i- j - J k-H 2 i- j - kx i 2 ,j,kExi 2 ,j,kH z i2 ,j 2 ,k， H z i2 ,j 2 ,k yn 1 iin =iiHy 2(i + 2 ,j,k %) -Hy 2 (i+ ,j,k-2A z上式中的第二项用平均值来替代 E： ”i +1, j,k)是因为离散方程中电场的时间取样是 整数n,磁场的时间取样是n+1/2,所以只能取n及n+1时电场的平均值。实际也证明这 个平均值使FD

15、TDJJ法具有数值稳定性。整理后，将EXG+j,k肝为未知数，其余作为迭代计算的已知数E?卡。+； ,j,k J=CA(m£(i +； ,j,k 计CB(m). 1111n11n11n11n11. Hz xAym=i,j,k+；（24）以上三式是电场的时间推进计算公式。2.5.2磁场时间推进差分格式节点（i,j,k用勺3个磁场分量分别用（i, j+,k + ）、（i +；, j, k + ； ）、（i +1 , j十,k位置上的Hx,Hy,Hz表示，同样，讨论式（4）中第一个公式，设观察点（x,y,z）为Hx的节点，即 在时刻t=nAt,对节点(i,j+1,k+2)的离散公式为：J1

16、 11211Hx 2 i, j 2,k 2 =CP m Hx 2 i, j 2,k 八七。+1,k+2 "E；(i,j,k + 2)E；(i,j +J,k + 1)-E”j+U-CQ(m )' (i +； ,j+； ,kHz 2 G +； ,j -1 ,k)_Hy 2 ( +；,j,k +；) Hy 2(i +； ,j,k yA zt一m =i ：, j,k；m m 1 m tCAm 一 ：t - 22； mCA m - ； m m - 1m ：t.：t 22 ； m.：tCB m =1: m；m m 一 1 . m ：tt 22 ； m同理，式(3)中其它两个公式的离散形式

17、为Ey+(, j +2,k "CA(mE；,j +2,k 计 CB(m>1111H：4(i, j +；,k + ； )-H：*(i, j +；,k ； ) H：'(i + ； , j + 2 ,kH：%i 一； , j +；,k)zAxm =i,j ：,kE；1 i,j,k -； =CAmE；i,j,k ； CBm-1111-II n q2 / ,1 . I . 1 II nd2 /1 . I . 1 I_jn 与1 1 I 1 1 II n42 /- -1 |上1(22)(23)Hy (i +2"卜2)Hy « -2，j，k、2)Hx (i, j

18、+2 卜十2)-H? (i, j - 2，k+ 2)(25)1 Im =i,j 2，kmm m1m m tCP m .，t - 2_ _2mmmm1mm.；ttCQ m =1 ； m:m1 m rt 22 ； m同理，式(4)中其它两个公式的离散形式为 11n o 11_n r 11Hy 2i 2,j,k 2 =CPmHy 2i 2，j,k2 E(i+；,j,k+1)-EXi +：,j,k) E；(i +1,j,k+6E；(i,j,k+：-CQ(m2222AzAxm=i+j,k+；(26)nJn 1H： 2i ；,j ；,k =CPmH：，i ；,j 晨k-E；(i+1,j +：kLE；G,j

19、 +：,k) E：(i +2，j +1,k)-EXi+2，j,k)1-CQ(m ) 2222AxAym = i +：，j +2*（27）以上三式是磁场的时间推进计算公式。时域推进计算框图（交叉半步逐步推进）在编程中，为了使电场和磁场有相同的数量级（为减小误差），可对H或E进行“归一化”处理，即：用H =Z°H取代H ,用 = £/2。取代E,式中Zo是自由空间波阻抗。计算结果再分别除以和乘以Z0即可。可以看出，这种离散方法电场和磁场在时间顺序上交替抽样，抽样间隔相差半个时间 步长，使麦克斯韦方程离散后成为显示差分方程，从而可以在时间上迭代求解，不需矩阵 求逆。给定初值后，可以逐步推进，求得以后各个时刻点的空间电磁分布。这是FDTD法的最大特点。2.6解的稳定性在FDTM,时间增量 组和空间增量Ax、Ay、Az之间不是相互独立的，它们的取值必须 满足一定的关系，以避免数值结果的不稳定，表现为随着时间步数的增加，计算结果发散 造成解不稳定的因素有多种： 误差因素：计算机在计算过程中，原始