第5章 CHAPTER5 时域有限差分法 5.1 时域有限差分法概述 时域有限差分法(FDTD)是由有限差分法发展出来的直接由麦克斯韦方程对电磁场进行 计算机模拟的数值分析方法,近年来时域有限差分法迅速地在电磁工程领域的各个方面得到 了广泛的应用和发展。本章将全面介绍时域有限差分法的有关理论、方法和应用。 5.1.1 时域有限差分法的特点 电磁场问题有3种基本的研究方法:理论分析、测量和计算机模拟。由于电磁场是以场 的形态存在的物质,具有独特的研究方法,很重要的特点就是要采取重叠的研究方法,也就是 只有理论分析、测量和计算机模拟的结果相互佐证,才可以认为是获得了正确可信的结论。时 域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。第4章已经详细地 探讨了有限差分法,有限差分法既能用于静电和时谐电磁场问题,也能用于时域电磁场问题的 数值分析。例如,第4章中介绍的抛物型和双曲型的偏微分方程中时间变量就是采用有限差 分法进行数值分析的。时域有限差分法的名称容易使人误解为只有时域有限差分法才具有处 理时域问题的能力,实际上,时域有限差分法的特点是直接由麦克斯韦方程组出发,在计算机 平台上对待求电磁学问题进行直接模拟,也可以认为是在计算机平台上构成的虚拟物理空间 上直接进行电磁过程的模拟试验。 近年来,许多学者对时域脉冲源的传播和响应进行了大量的研究,谈到时域电磁学问题, 首先要关注的应该是描述物体在瞬态电磁源作用下普遍适用的理论。例如,奇点展开方法 (SingularityExpansionMethod,SEM)就是关于如何利用复频域上物体对电磁波响应的极点 来描述物体在瞬态电磁源作用下响应的普遍适用的理论。该方法是对电磁脉冲响应进行数值 模拟的理论基础,它证明了瞬态电磁源作用下响应的极点或奇异点只发生在引起物体自然谐 振的复频率点上。自然谐振点只与物体的几何形状和几何尺寸有关,而与外界的激励源无关。 所以,时域有限差分法只是处理时域电磁学问题的一种方法而不是唯一的方法。此外,尽管时 域有限差分法是直接由麦克斯韦方程组出发的,但是决不能取代电磁学理论分析,时域有限差 分法的结果必须能够禁得起理论的考验,当两者的结论发生矛盾时,一般都可以归结为数值处 理过程不严密或发生错误。 物体的电特性基本是时谐的,物体的电参数都是频率的函数。世界上真实的时域电信号 总是具有有限能量的(时谐信号除外),理论上具有几乎所有频率成分,但是对实际的电磁学问 163 题来说,只有有限的频带内的频率成分在起主要作用。所以,实际上不需要考虑真正的具有无 限大频谱成分的时域电信号,只需要考虑具有一定带宽的能量信号,更接近于具有群速度和有 限带宽的时谐波群的情形。 认为时域有限差分法是时域方法而其他方法都是频域方法是不正确的,也就是说使用任 何一种计算电磁学方法都需要考虑所处理问题的尺寸和问题关注的电磁波的波长,区别只是 由频域方程出发的方法需要配合快速傅里叶变换求时域问题的解。计算电磁学方法本身并没 有明确的频率范围的限制,对具体的电磁学问题来说,采用什么计算电磁学方法要根据物体的 电尺寸来决定。物体的电尺寸决定了物体的电谐振频率,所以通常以研究的物体尺寸所决定 的谐振频率为参考点,将计算电磁学方法分为高频、中频和低频方法。以下是使用计算电磁学 方法时要考虑的因素。 首先,计算机内存和计算需要的CPU 时间限制了计算电磁学方法的应用范围。采用时 域有限差分法、矩量法和有限元法时所产生的矩阵大小是限制应用的关键因素,以时域有限差 分法为例,目前可以达到104×104×104 个单元的尺度。其次,各种方法的数值建模的处理方 法也限制了应用范围。最后,在物理问题进行数值化时,基本的要求是要保持物理问题和物理 规律在数值化的过程中不发生明显改变,所以量化时产生的数值色散等问题,也决定了不同方 法的不同应用范围。 通常矩量法和有限元法只能分析谐振点以下的问题,可以称为低频方法;时域有限差分法 用来分析系统第一谐振点频率上下谐振波幅度4个量级的范围,可以称为中频方法;几何射线 法用来分析远高于谐振点的电磁学问题,可以称为高频方法。 FDTD 法有如下特点。 (1)适用于分析系统谐振点附近很宽的频带响应。 (2)可以分析任意三维形状的问题。 (3)适用于研究理想导体、实际金属和绝缘物体等各类物体在电磁波作用下的效应。 (4)适用于处理具有频谱依赖性的媒质参量,如损耗介质、磁介质、非寻常物质(如各向异 性媒质、铁氧体)等的电磁学问题。 (5)适用于分析任意类型的响应,包括远场和近场,如散射场、天线方向图、雷达散射截面 (RCS )、表面波、电流、功率密度、穿透和内耦合等。 (6)适用于分析雷电、EMP 、HPM 、雷达、激光器等激励源。 (7)适用于分析多种多样的系统,如烟雾、屏蔽或防护罩、飞机、人体、卫星、探测等。 采用混合方法(hybridtechniques), 例如GTD 或DO+FDTD 可以分析低于谐振频率直 到甚高频范围的电磁问题。总之,FDTD 可分析系统的第一谐振频率以上谐振波幅度为4个 量级的范围的问题,也就是能达到由低到高频谐波幅度达到6个数量级的范围,按功率约合 120dB,按场强约合60dB 的范围。 设工作站上能处理约100 万(106)个单元,对三维问题来讲,具有约100×100×100 个单 元大小。若按标准的情况,每波长取10 个空间步长,这个工作站能处理约10 个波长立方空间 中的电磁问题。现在的超级计算机每小时可处理约107 个单元,约达到46 个波长的空间,也 就是说,不同硬件配置的计算机所能创造的数值空间大小是有不同限度的,解决问题的能力也 不同。 总之,FDTD 具有能处理宽频带、多种模拟源、各种形状作用物和复杂环境的电磁学问题 的能力,具有采用计算机类型宽、响应量级宽等方面的优势;具有计算效率高并且可以直接得 1 64 到宽频带结果的优势。因此,FDTD特别适于处理薄板和细导线天线问题,只要计算机能够处 理足够多的FDTD单元,计算的精度就可达到任意要求。此外,FDTD程序可以规范化、商品 化,易于推广,与图形工具结合可以很方便地用来分析电磁学问题,因此它实际上是一种集实 验方法、计算方法和分析方法于一身的方法。 在如下常见的6类问题中,除第(1)类问题之外,都可以采用FDTD方法处理。 (1)电源———电力、设备等问题。 (2)传输线、波导等传输问题。 (3)天线的接收、检测和辐射。 (4)耦合、屏蔽和透入效应。 (5)散射和逆散射问题。 (6)开关、过渡过程等非线性问题。 FDTD最适宜分析的是瞬态响应问题,特别是具有复杂几何形状和复杂环境的情形,例如 埋地天线、介质覆盖天线等。MoM 在分析高频响应问题时,往往误差过大,特别是用MoM 分 析封闭金属体内接近谐振点的问题时,会产生很大的误差,此时采用FDTD就很合适。FDTD 很难分析低频响应,例如电力线的传输问题如果采用FDTD就要求很多时间步,甚至需要10 亿个时间步,此时采用MoM 应该是很好的选择。 5.1.2 电磁场旋度方程 为了介绍FDTD的基本原理,首先分析线性介质中麦克斯韦方程的特点。在线性介质 中,麦克斯韦方程可以表达为式(5.1.1)的形式。 Δ×E =-.B .t Δ×H =.D .t +J Δ·D =ρ Δ·B =0 ì . í .... .... (5.1.1) 式中 B =μH D =εE 通常在时间起点,场和源都置0,此时,两个散度方程可以包括在两个旋度方程和初始的边界 条件之中,是冗余的,FDTD公式只需要麦克斯韦方程组的旋度方程就足够了,即 .H .t =-1μ ( Δ×E)-σ* μ H .E .t =σ ε ( Δ×H )-σ εE 此处,J=σE,也应包含磁损耗介质,故采用σ* 表示磁导率。FDTD中的变量采用电场强度E 和磁场强度H ,不是磁感应强度B 和电位移矢量D。下面的推导可以说明散度方程已经包含 在旋度方程之中,取 Δ· Δ×E =-.B .t . è . . . ÷ 因为Δ · Δ ×A≡0,所以有Δ ·B=常数。 1 65 取 Δ· Δ×H =.D . .t +J è . . . ÷ 得 0=-.( Δ·D) .t + Δ·J 由电荷守恒定律,有 Δ·J +.ρ .t =0 于是 .( Δ·D) .t -.ρ .t =0 因而有 . .t[( Δ·D)-ρ]=0 必有 Δ ·D-ρ=常数 因为在FDTD计算时,开始时场和源都置0,所以 Δ·B t=0=0, Δ·D -ρ t=0=0 由这两个条件就可以定出这两个常数为0,则有 Δ·B =0, Δ·D -ρ=0 于是,由旋度方程推出了散度方程,也就证明了散度方程已经包含在旋度方程之中了。虽然两 个散度方程在FDTD中不出现,但是在实际应用时,仍然有用,可以用它们来验证最后结果的 正确性。 5.1.3 分裂场形式 因为麦克斯韦方程组是线性方程组,电磁场可以分裂为入射场和反射场之和的形式,即 E =Etotal ≡Einc+Escat H =Htotal ≡Hinc+Hscat { (5.1.2) 这样分裂的好处是:入射场可以在整个问题空间中用解析法求解,而只有散射场需要进 行数值分析,只有散射场需要在研究的问题空间之外被边界条件吸收。这一点很重要,因为散 射场要比全场更容易设置吸收边界条件。 入射场在自由空间单独地、独立地满足麦克斯韦方程组,即 Δ×Einc=-μ0.Hinc .t Δ×Hinc=ε0.Einc .t ì . í .. . .. (5.1.3) 由麦克斯韦方程组的旋度方程,一般有 Δ× (Einc+Escat)=-μ.(Hinc+Hscat) .t -σ* (Hinc+Hscat) Δ× (Hinc+Hscat)=ε.(Einc+Escat) .t +σ(Einc+Escat) ì . í .. . .. 1 66 将自由空间条件下的入射场方程代入,有 Δ×Escat=-μ.Hscat .t -σ*Hscat- (μ -μ0).Hinc é .t +σ*Hinc . êê ù . úú Δ×Hscat=ε.Escat .t +σEscat- (ε-ε0).Einc é .t +σEinc . êê ù . úú ì . í .. . .. (5.1.4) 在散射区外趋向无限远处的电磁过程可以看成与自由空间的情形一样,也就是说入射场与散 射场的和满足自由空间中的麦克斯韦方程组,即 Δ×Etotal=-μ0.Htotal .t Δ×Htotal=ε0.Etotal .t ì . í .. . .. 所以,在远场有 Δ× (Einc+Escat)=-μ0.(Hinc+Hscat) .t Δ× (Hinc+Hscat)=ε0.(Einc+Escat) .t ì . í .. . .. 代入式(5.1.3),可得到在自由空间中,有 Δ×Escat=-μ0.Hscat .t Δ×Hscat=ε0.Escat .t ì . í .. . .. (5.1.5) 实际上,在式(5.1.4)中,令μ=μ0,ε=ε0,σ→0和σ* →0,替换后也可得到式(5.1.5)。所以,入 射场是自由空间的电磁学问题,只有散射场必须满足式(5.1.4),必须用数值方法求解。式(5.1.4) 写为 .Hscat .t =-σ* μ Hscat-σ* μ Hinc- (μ -μ0) μ .Hinc .t -1μ ( Δ×Escat) .Escat .t =-σ εEscat-σ εEinc- (ε-ε0) ε .Einc .t -1ε ( Δ×Hscat) ì . í ... .. . (5.1.6) 5.1.4 理想导体的时域有限差分法公式 在散射体外,散射场满足自由空间条件,即 σ* =σ=0 μ =μ0 ε=ε0 则式(5.1.6)变为 .Hscat .t =- 1 μ0( Δ×Escat) .Escat .t =1 ε0( Δ×Hscat) ì . í .. . .. (5.1.7) 由式(5.1.6),在导体中,有 ε σ .Escat .t =-Escat-Einc- (ε-ε0) σ .Einc .t +1σ ( Δ×Hscat) 1 67 因理想导体条件为σ=∞,可简化为 Escat=-Einc (5.1.8) 理想导体为散射体时,电磁场可以用式(5.1.7)描写。因此,如果问题中只有自由空间和 理想导体,则采用FDTD时,只需式(5.1.7)和式(5.1.8)就足够了。 把散射场表示为分量的方程,有 TE波 .Escat x .t =1 ε0 .Hscat z .y -.Hscat y .z . è . . . ÷ .Escat y .t =1 ε0 .Hscat x .z -.Hscat z .x . è . . . ÷ .Hscat z .t = 1 μ0 .Escat x .y -.Escat y .x . è . . . ÷ ì . í .... .... , TM 波 .Hscat x .t = 1 μ0 .Escat y .z -.Escat z .y . è . . . ÷ .Hscat y .t = 1 μ0 .Escat z .x -.Escat x .z . è . . . ÷ .Escat z .t =1 ε0 .Hscat y .x -.Hscat x .y . è . . . ÷ ì . í .... .... (5.1.9) 用差分式代替微分式,有 Escat,n x -Escat,n-1 x Δt =1 ε0 ΔHscat,n-12 z Δy -ΔHscat,n-12 y Δz . è . . . ÷ Hscat,n+12 y -Hscat,n-12 y Δt = 1 μ0 ΔEscat,n z Δx -ΔEscat,n x Δz . è . . . ÷ ì . í ... .. . (5.1.10) 5.1.5 损耗媒质的情况 将分裂场的方程代入式(5.1.6),有 E =Etotal ≡Einc+Escat H =Htotal ≡Hinc+Hscat { 即 .Hscat .t =-σ* μ Hscat-σ* μ Hinc- (μ -μ0) μ .Hinc .t -1μ ( Δ×Escat) .Escat .t =-σ εEscat-σ εEinc- (ε-ε0) ε .Einc .t -1ε ( Δ×Hscat) ì . í ... .. . (5.1.11) 有损耗媒质中的电磁场的基本关系式(5.1.11)可以简化为 ε.Es .t +σEs =-σEi - (ε-ε0).Ei .t - Δ×Hs (5.1.12) 式中,采用Es 表示Escat,Ei 表示Einc,其他变量也类似地简化书写。 5.2 时域有限差分法基础 5.2.1 使用时域有限差分法的影响因素 在使用FDTD方法时必须处理单元尺寸、时间步长、入射场、散射体结构、场强计算、吸收 边界条件及资源需求等问题。在FDTD 计算中,确定单元尺寸是很关键的步骤,单元尺寸必 须足够小,应该使最高频率的结果有足够的精度,当然所确定的单元尺寸不能太小,不能超出 计算机的资源允许的计算能力范围,必须是可实现的。媒质参数也影响空间步长:ε、μ 和σ 大,则波长短,给定频率下的单元尺寸就应当更小。 1 68 空间步长确定之后,就可由数值计算稳定条件决定时间步长。当采用FDTD 时,入射场 应当是能够解析表达的,入射场的频率特性对时间步长影响很大。许多问题允许选择多种形 式的入射场信号,通常只要可能就将入射场波形设为高斯脉冲,因为高斯脉冲的带宽最窄。当 媒质特性与频率有关时,选平滑余弦脉冲(smoothedcosinepulse)源就较好。入射场源的处 理是实现FDTD的另一个关键问题,当然也应该设定恰当的入射场源馈入的机制。 FDTD能达到0.1dB的精度和120dB的动态范围,实际计算效果取决于单元尺寸、频率和 物体形状,决定于采用的FDTD方程和系数的取值,决定于如何求取关联变量的数值。 恰当地设定吸收边界条件(absorbingboundaries,absorbingboundaryconditions)是实现 FDTD算法的另一关键问题,一般情况都采用MUR一阶和二阶吸收边界,当然采用其他种类 的吸收边界可以获得更好的吸收效果,但相应地也增加了程序的复杂度和计算量。 时间步数应该足够多,要能显示出作用场的特性,特别是要使数值计算的结果能够给出谐 振特性。在用FDTD方法讨论时谐场问题时一定要有足够的计算步数,直到得到稳定的周期 解为止。 当确定了问题中Yee单元的个数和时间步数之后,就可以估计所需要的CPU 时间和占 用内存(RAM)量、所需硬盘存储空间和总花费等。 5.2.2 Yee 单元网格空间中电磁场的量化关系 在推导FDTD差分格式时,采用中心差商代替微商并且用正六面体网格进行空间切分, 产生的量化空间(quantizespace)具有如下关系: x =IΔx, y =JΔy, z =KΔz, t=nΔt 此时,用I、J、K 就可以表示网格空间的坐标。取均匀六面体网格和中心差商并考虑电磁场 分量之间的方向和旋度关系,就导致了图5.1中的Yee单元网格。 图5.1 Yee单元网格空间中的电磁场 图5.1中,En z(I,J,K )表示电场的Z 分量在x=IΔx,y=JΔy,z= K +12 . è . . . ÷ Δz,t=nΔt 的值,其他分量也可类似理解。Yee单元的建立是采用FDTD算法由麦克斯韦方程出发进行 电磁过程模拟的关键,该单元把数学关系、物理含义和物理规律巧妙地结合在一个差分单元 中,Yee单元本身就是一套麦克斯韦方程,因此应该从有限差分和电磁定律相结合的角度进行 1 69 理解。图5.1和图5.2分别从两个角度画出了Yee单元,其特点如下: 图5.2 Yee单元网格及场关系 (1)电场与磁场分量在空间交叉放置,相互垂直。 (2)每个坐标平面上电场分量的四周由磁场分量环绕,磁场分量的四周由电场分量环绕。 (3)每个场分量自身相距一个空间步长,电场与磁场相距半个空间步长。 (4)电场取n 时刻、磁场取n+1/2时刻的值。 (5)电场的n+1时刻的值由n 时刻的值得到,磁场的n+1/2时刻的值由n-1/2时刻的 值得到,电场的n+1时刻的旋度取n 时刻电场的值,磁场的旋度取n-1/2时刻磁场的值(对 应H 为n+1/2时刻的旋度)。对应于电场旋度产生磁场、磁场旋度产生电场的定律,不能违 反因果率。例如,电场的n 时刻的旋度产生的应该是下一时间步的磁场,也就是n+1时刻电 场所对应的磁场,即(n+1)+1/2时刻的磁场。 (6)由于时间为中心差商,Yee单元网格内的空间亦为中心差商。由于均匀媒质中电磁 波的空间变量和时间变量完全对称,所以在Yee单元网格内,电场和磁场的空间位置也应该 差半个空间步,而且在3个空间坐标方向上的时间步必须相等。 (7)Yee单元网格内的媒质只取一种。单元递推时,如果需要用到相关单元的媒质,要取 单元网格自身的和最邻近单元的媒质参数值。简言之,Yee单元网格内完全用均匀媒质和平 面电磁波近似,Yee单元网格内电场与磁场是不变的,只保留着定律的关系。 (8)因为Yee单元本身就是一组数值化和几何化了的麦克斯韦方程,所以Yee网格单元 在数值建模和计算中是不可拆卸的,在Yee单元外部的单元布局与有限差分法是类似的,要 根据问题的需要决定,但是不管什么情况下,一定要保证Yee单元内部关系的完整性,否则就 破坏了麦克斯韦方程的关系,就会产生错误。请读者考虑,如果面对柱面波或球面波应该如何 处理Yee单元网格? 采用立方体的Yee单元网格与采用柱面或球面的Yee单元网格有什么 区别? 5.2.3 决定单元的空间尺寸 在应用FDTD时,选择单元的空间尺寸是很重要的。单元的尺寸一定要小于最短波长, 为了保证计算的精度,人们认为选取单元尺寸为λ/10是一种规则,实际上,单元尺寸可以根据 1 70 具体情况选择,并不是一成不变的。以下是一些选取单元尺寸的考虑。 (1)通常选每波长10个单元,即选单元尺寸为λ/10,或更小。在需要更高精度的场合,例 如决定雷达散射截面的情形需要选取单元尺寸为λ/20或更小。 (2)原则上,只要满足一个波长中有4个空间单元就能得到合理的结果,这是由奈奎斯特 (Nyquist)抽样准则决定的,在信号处理中,奈奎斯特抽样准则用于时间抽样,而现在是将奈奎 斯特抽样准则用于空间抽样,其原理完全一样。奈奎斯特抽样准则要求λ=2Δx,即要求在一 个波长(对应时间为周期)中至少要有两个抽样,如果单元尺寸比要求的抽样间隔小得多,模拟 将与真实电磁波的情况非常接近,就一定能得到合理的结果。但实际上通常要选更小的单元, 因为每个时间步网格都同时是电场和磁场在Yee单元网格空间的抽样,所以一个波长范围最 好要抽样4次以上。 (3)由于网格抽样并不精确,也无法预先精确得到研究问题中的最小波长,因而就会存在 网格发散误差(griddispersionerror)。由于FDTD方法本质上带来的近似,不同频率的波在 网格中传播时会有不同的速度(物理上是同一速度),差异的程度取决于传播方向与网格方向 的关系。因此,为了使误差控制在可接受的水平上,就需要更小的单元尺寸,显然每波长抽样 4次是远远不够的。 (4)与单元尺寸有关的另一因素是问题的几何形状和电尺寸。通常只要取单元尺寸为 0.1λ 或更小,就总可以满足问题在几何形状上的要求。但是,一些特殊的几何形状会要求更 小的网格尺寸,就应该在根据研究问题的频率确定单元网格尺寸的同时照顾到几何形状的要 求。典型的情形是设计线天线时,虽然天线的直径或厚度为λ/10到λ/20或更小,但是却直接 影响天线阻抗。又例如,对圆形物体建模,采用矩形逼近边界时,阶梯效应会引起显著的误差。 在这些情况下就需要采用更小的网格或采用特殊方法进行处理,例如采用子单元模型或采用 比Yee单元能更好模拟实际情况的单元模型。 单元尺寸决定后,物体需要多少个单元、物体周围空间需要多少个单元也就确定了, FDTD的Yee单元网格空间的大小也就确定了,一般的三维问题总要有数百、数千、数百万个 单元。据此,就可以结合使用的计算机(如PC、工作站和并行机)决定所需要的计算时间和 花费。 5.2.4 离散化的麦克斯韦方程 1.无源无耗媒质中的FDTD 格式 在无源无耗媒质中,将式(5.1.1)中的磁场旋度方程按照标量形式展开,得 ε.Ex .t =.Hz .y -.Hy .z ε.Ey .t =.Hx .z -.Hz .x ε.Ez .t =.Hy .x -.Hx .y (5.2.1) 采用图5.2所示的Yee单元网格对式(5.2.1)表示的麦克斯韦微分方程进行差分。设Δx、 Δy 和Δz 分别代表在x、y 和z 坐标方向的空间步长,在第n 步,采用二阶中心差分形式,例如 对于Ez ,则在x 方向及时间t 上的差分格式分别为 1 71 .En z (i,j,k) .x = Enz i+12 ,j,k . è . . . ÷ -Enz i-12 ,j,k . è . . . ÷ Δx +O(Δx2) .En z(i,j,k) .t = En+12 z (i,j,k)-En-12 z (i,j,k) Δt +O(Δt2) 其中,i、j、k 分别表示计算空间沿x、y、z 坐标方向的步数。经过以上形式的差分近似,式(5.2.1) 有如下形式: En+1 x i+12 ,j,k . è . . . ÷ =En x i+12 ,j,k . è . . . ÷ + Δt ε Hn+12 z i+12 ,j+12 ,k . è . . . ÷ -Hn+12 z i+12 ,j-12 ,k . è . . . ÷ Δy é . êêêêê - Hn+12 y (i+12 ,j,k +12 )-Hn+12 y (i+12 ,j,k -12 ) Δz ù . úúúúú En+1 y i,j+12 ,k . è . . . ÷ =En y i,j+12 ,k . è . . . ÷ + Δt ε Hn+12 x i,j+12 ,k +12 . è . . . ÷ -Hn+12 x i,j+12,k -12 . è . . . ÷ Δz - é . êêêêê Hn+12 z i+12 ,j+12 ,k . è . . . ÷ -Hn+12 z i-12 ,j+12 ,k . è .. . ÷ Δx ù . úúúúú En+1 z i,j,k +12 . è . . . ÷ =Enz i,j,k +12 . è . . . ÷ + Δt ε Hn+12 y i+12,j,k +12 . è . . . ÷ -Hn+12 y i-12 ,j,k +12 . è . . . ÷ Δx - é . êêêêê Hn+12 x i,j+12 ,k +12 . è . . . ÷ -Hn+12 x i,j-12 ,k +12 . è . . . ÷ Δy ù . úúúúú 同理,将麦克斯韦方程中电场旋度方程的微分形式进行差分,得 Hn+12 x i,j+12 ,k +12 . è .. . ÷ =Hn-12 x i,j+12 ,k +12 . è . . . ÷ + Δt μ En y i,j+12 ,k +1 . è . . . ÷ -En y i,j+12 ,k . è . . . ÷ Δz - é . êêêêê