应用FDTD方法解决电磁辐射问题

应用FDTD方法解决电磁辐射问题自电磁场基本方程以来，电磁场理论和应用的发展已经有一百多年的历史。

目前，电磁波的研究已深入到各个领域，应用十分广泛，例如无线电波传波，光纤通信和移动通信，雷达技术，微波，天线，电磁成像，地下电磁探测，电磁兼容等等。

在各类复杂系统中的电磁问题，主要依靠各种电磁场数值计算方法加以解决。

随着电子计算机处理能力和存储容量的巨大发展，更促进了这些计算方法在实际问题中的应用。

目前在电磁场领域应用的数值算法也是种类繁多，各有其优缺点，常用的电磁场计算方法大致有：FDTD Finite difference time domain （时域有限差分法）TLM Transmission line method （传输线法）FEM Finite element method （有限元法）BEM Boundary element method （边界元法）MoM Method of moments （矩量法） 其中时域有限差分法(FDTD)理论经过30多年的发展和完善，已经成为时域电磁场数值计算的主要方法之一，并广泛应用各类实际工程电磁场中。

一、 FDTD 法简介时域有限差分法以差分原理为基础，直接从概括电磁场普遍规律的麦克斯韦旋度方程出发，将其转换为差分方程组，在一定体积内和一段时间上对连续电磁场的数据采样。

因此，它是以电磁场问题的最原始、最本质、最完备的数值模拟。

以它为基础制作的计算程序，对广泛的电磁场问题具有通用性，因此得到了广泛的应用。

1. Yee 差分算法基本原理考虑空间一个无源区域，其煤质参数不随时间变化且各向同性，由Maxwell 方程组中的两个旋度方程在直角坐标系中可导出六个耦合公式：1(1.1)1(1.2)H E H t E H E t ρμμσεε∂=-∇⨯-∂∂=∇⨯-∂ ⇒1()(1.3)1()(1.4)1()(1.5)1()(1.6)1()(1.7)1()(1.8)y x z x y x z y y x z z y x z x y x z y y x z z E H E H t z y H E E H t x z E E H H t y x H E H E t y z E H H E t z x H H E E t x y ρμρμρμσεσεσε∂⎧∂∂=--⎪∂∂∂⎪⎪∂∂∂=--⎪∂∂∂⎪⎪∂∂∂⎪=--∂∂∂⎪⎨∂∂∂⎪=--⎪∂∂∂⎪∂∂⎪∂=--⎪∂∂∂⎪∂∂∂⎪=--⎪∂∂∂⎩其中ε为介电常数（F/m ）；μ为磁导率（H/m ）；σ为电导率（S/m ）；ρ为磁阻率（/m Ω）。

按照Yee 差分算法，首先在空间建立矩形差分网格，网格节点的空间坐标与一组相应的整数标号一一对应： (,,)(,,)i j k i x j y k z =∆∆∆而该点的任意函数F(x,y,z,t)在时刻n t∆的值可以表示为(,,)(,,,)(1.9)n F i j k F i x j y k z n t =∆∆∆∆ ,,x y z∆∆∆分别为矩形网格沿x ，y ，z 方向的空间步长t ∆是时间步长。

时域差分法的实质就是将计算问题在时间和空间进行离散: 用中心差分取二阶精度：对空间离散：211(,,)(,,)(,,)22()(1.10)n n n F i j k F i j k F i j k O x x x +--∂=+∆∂∆ 对时间离散：11222(,,)(,,)(,,)()(1.11)n n n F i j k F i j k Fi j k O t t t +-∂-=+∆∂∆ 为了获得（）式的精度，并满足（）至（），可将空间任意矩形网格上的E 和H 的六个分量如图1-1所示放置。

图1-1 Yee 网格单元及电磁场空间离散点关系 每个磁场分量由四个电场分量环绕着；每个电场分量由四个磁场分量所环绕。

为获得（）的精度，可将E 和H 在时间上相差半个步长交替计算。

按照这些原则，可将式（）~（）化为差分方程组。

以（）式为例，差分后得到（）1111122221(,,)2112(,,)1112(,,)(,,)11122(,,)(,,)(,,)22211112(,,)2(,,)221111111(,,)(,,)(,,)(222222n n x x n n n n z z y y i j k ti j k t E i j k E i j k i j k t i j k i j k t i j k i j k H i j k H i j k H i j k H i y σεσεσεε++++++∆-+∆+=⋅++⋅+∆++∆++++++-+-++-+⋅+∆1,,)22(1.12)j k z ⎡⎤-⎢⎥⎢⎥∆⎢⎥⎢⎥⎣⎦其余的也可以写出，每个网格点上的个场分量的新值依赖于该点在前一时间步长时刻的值即该点周围临近点上另一场量在早半个时间步长时的值。

因此任意时刻可一次算出一个点，并行算法可计算出多个点。

通过这些运算可以交替算出电场磁场在各个时间步的值。

具体参考[1][2]，这里不多叙述。

2. 数值稳定性条件与网格划分数值计算中一个比较重要的问题就是数值解的稳定性问题，这在FDTD 方程中取决于网格的空间步长和时间不长的关系。

数值稳定性条件以上所建立的麦克斯韦旋度方程的有限差分算式是显式的，时间步长t∆，空间步长,,x y z∆∆∆必须满足一定关系，否则会使得数值表现不稳定。

这种不稳定表现为在解显示差分方程时，随着时间步数的增加，计算结果也无限制的增加。

Taflaove 等对Yee差分格式稳定性进行讨论，导出了对时间步长的限制性条件：(1.21)t∆≤当x y z s∆=∆=∆=∆时，数值稳定条件（12）可简化为(1.22)t s∆≤=一维条件下，这一稳定性条件就简化为：(1.23)stv∆∆≤上式可以看出，即要求时间步长不能大于电磁波传播一个空间步长所需的时间，否则，就破坏了电磁波传播的因果关系[12]。

如果计算空间中的媒质是不均匀的，那么稳定条件对不同的煤质区域是不同的，我们只需取最大的V 满足的条件，在其他区域中也就得到满足了。

此时有： max (1.24)t n ν∆≤ 其中n 为空间维度。

网格划分由数值稳定性条件可知，网格划分需要满足一定的条件，在有足够条件的情况下，将网格划分的越细，计算的结果就会越精确。

目前有一款交互式软件网格划分软件GEOM ，根据设计需要，可自动的用FDTD 的方法分析几何结构，产生离散模型。

但是计算机计算能力毕竟是有限的，有时候提高一定比例的网格大小，并不一定能够获得更多的精度，同时计算的效率却下降了很多。

而且对于一些复杂的模型结构，由于介质层的物理尺寸很小，较大的网格无法覆盖这样的结构，就导致了划分后介质层的缺损。

因此需要在满足稳定性条件的前提下对划分的大小进行优化。

最好的办法就是采用非均匀网格划分的办法。

合理的利用了计算区域内不同的计算需要，对于一些关键位置的地方，如不同介质的交接处，或者是结构细微处，或者是馈源的地方，对网格细化，而在介质均一的位置，则放大网格的大小，使得总的网格数量得到降低。

对于非均匀划分方案的确定，也并不是盲目的过程。

可以首先利用较大的网格大小进行一下预模拟，对模拟的结果进行分析，判断电场(或磁场等其他场分布)分布的规律及形态，对于关键计算的部位，再采用较小的网格划分。

其次，可以利用遗传算法等优化算法，对网格划分进行局部最优化的搜索，提高网格划分的质量。

下图是对手机天线几种网格划分示意图。

图均匀网格划分图非均匀网格划分图含有亚网格的网格划分除了非均匀网格划分的方法，另外一种在这个基础上的改进网格划分方法就是亚网格划分的方法[14]。

亚网格划分在原有网格划分的基础上，在计算的关键点，或者是重点关注的部分，或者是介质层物理尺寸很薄的地方，将网格进一步划小，而其它部分的网格则保持不变。

这样做的好处在于可以在不影响其它部分的情况下，提高重点计算区域的计算精度。

同时网格的数量不会因此激增。

3.吸收边界条件当我们考虑像辐射和散射这种开放向题时，原则上应该在无限空间中计算电磁场，但是，由于计算机存储量的限制不可能考虑很大的空间,，因而，一个基本问题是如何处理网格的截断向题。

一种有效的方法是在截断处设置一种吸收边界条件，如果它能吸收全部向外传播的波，就等效于考虑了一个无限大的空间。

目前构造吸收边界条件的思路主要有两种：一种是从外形波方程出发构造的透射边界条件，如Mur边界条件，这种类型的透射边界条件具有构造简单，内存需求小等特点。

另一种是在边界上引入吸收材料，电磁波在无反射地进入吸收材料后被衰减掉，如PML。

这种方法构造复杂，内存需求量较大，但在很大的入射角度的吸收条件较好。

Mur给出了适合于在时域有限差分应用的吸收边界条件的一阶和二阶的近似形式，具体推导过程可见[16][17]。

应该注意的是Mur的局限性：（1）边界面上任意点的插值都是在其领域的三维空间上进行的，不允许靠近边界的区域出现介质的不连续；（2）在垂直入射的情况下的计算精度仅相当于一阶吸收边界条件。

完全匹配层PML是通过在截断边界外放置虚拟的吸收介质来系数外形波的，PML边界由硬边界构成，在计算PML边界时，就可以不涉及边界外的未知量，大幅度提高了计算精度，详细过程见[18]与Mur边界条件相比，PML边界具有更高的计算精度而且任意方向入射时均无反射。

在二维自由空间，PML边界吸收的反射系数可低于-70dB，比其他各种边界条件改善约40dB。

这种人工涉及的完全匹配层由有耗导电、导磁煤质组成，可吸收人以入射角、任何频率、任意偏振态的入射电磁波，实用性更强，它内存大，尽管有的文献[19]提出的各向异性媒质的PML 吸收边界条件所需计算机内存量减少了一半，但是其实际计算量并没有减少多少，在不连续介质接触边界面和三维区域的6个面不全为吸收边界时也难以应付，具体算法可见[19]。

目前各种吸收边界条件都存在这一些缺点，吸收边界条件反射过大，所适应的入射角的范围不够宽广，提高计算精度等都应该是我们努力的方向。

二、电磁辐射问题电磁辐射是由空间共同传送的电能量和磁能量所组成，而该能量是由电荷移动所产生的。

时变的电荷和电流是激发电磁波的源。

为了有效使电磁波能量按所要求的方向辐射出去，时变的电荷和电流必须按某种特殊的方式分布，天线就是设计成按规定的方式有效的辐射电磁波能量的装置。

其实质就是一种导行波与自由空间的转换器件或者换能器，用于发送和接受自由空间的电磁波。

举例说，正在发射讯号的射频天线所发出的移动电荷，便会产生电磁能量，电磁包括形形色色的电磁辐射，从低频电磁辐射到高频电磁辐射，两者之间还有无线电波、微波、红外线和紫外光等。

任何形式和形状的天线都可以看成是有无限多个载有交流信号的基本小线元组成，这些基本元通常被称为电偶极子，电偶极子产生的电磁场分析计算是线性天线工程的基础。