有限差分法教材课件.ppt

《有限差分法教材课件.ppt》由会员分享，可在线阅读，更多相关《有限差分法教材课件.ppt（56页珍藏版）》请在三一办公上搜索。

1、第 4 章 有限差分法,本章基于差分原理阐述了在电磁场数值计算方法中应用最早的有限差分法，并以正方形网格划分的离散模式为主体，重点讨论了静态场中方法应用的全过程，并介绍了时变电磁场中直接将麦克斯韦方程组中的旋度方程转化为差分方程的时域有限差分法。,4.1 概述,在电磁场数值计算方法中，有限差分法(Finite Difference Method，简称FDM)是应用最早的一种方法。有限差分法以其概念清晰，方法简单、直观等特点，在电磁场数值分析领域内得到了广泛的应用。现阶段各种电磁场数值计算方法发展很快，尤其是在有限差分法与变分法相结合的基础上形成的有限元法日益得到广泛的应用，但有限差分法以其固有

2、的特点仍然是一种不容忽视的数值计算方法。例如，面向高频电磁场的传输、辐射、散射和透入等工程问题的需求，基于麦克斯韦方程组中旋度方程直接转化为差分方程的时域有限差分法(Finite Difference Time Domain Method，简称FDTD)即从传统的有限差分法中脱颖而出，成为在上述一系列工程问题中广泛应用的数值计算方法。,为求解由偏微分方程定解问题所构造的数学模型，有限差分法的基本思想是利用网格剖分将定解区域（场域）离散化为网格离散节点的集合，然后，基于差分原理的应用，以各离散点上函数的差商来近似替代该点的偏导数，这样，待求的偏微分方程定解问题可转化为相应的差分方程组（代数方程组

3、）问题，解出各离散点上的待求函数值，即为所求定解问题的离散解，若再应用插值方法，便可从离散解得到定解问题在整个场域上的近似解。,对于包括电磁场在内的各种物理场，应用有限差分法进行数值计算的步骤通常是：,1）采用一定的网格剖分方式离散化场域；2）基于差分原理的应用，对场域内偏微分方程以及定解条件进行差分离散化处理(一般把这一步骤称为构造差分格式)；3）由所建立的差分格式(即与原定解问题对应的离散数学模型代数方程组)，选用合适的代数方程组的解法，编制计算程序，算出待求的离散解。,有限差分法有上述大致固定的处理和计算模式，具有一定的通用性。,4.2 差分与差商,有限差分法是以差分原理为基础的一种数值

4、计算法。它用离散的函数值所构成的差商来近似逼近相应的偏导数， 而所谓差商则是基于差分应用的数值微分表达式。 设一函数 f(x)， 其自变量 x 得到一个很小的增量x = h， 则函数 f(x)的增量,称为函数 f（x）的一阶差分。显然，只要增量 h 很小， 差分f与微分 df之间的差异将很小 。,一阶差分仍是自变量 x 的函数，相类似地按式（4-1）计算一阶差分的差分， 就得到2f(x)，称之为原始函数 f(x)的二阶差分。 同样， 当 h 很小时， 二阶差分2f(x)逼近于二阶微分d2f。依同理，可以定义更高阶的差分。,即是无限小的微分 除以无限小的微分 的商，应用差分，显然，它可近似地表达

5、为,即有限小的差分f(x)除以有限小的差分x 的商，称为差商。同理，一阶导数 还可近似表达为,一阶导数,式（4-2）、 式（4-3）和式（4-4）分别称为一阶向前、 向后和中心差商。 如图 4-1 所示， 对应于点 P 的一阶向前、 向后和中心差商，在几何意义上可分别表征为弧线 PB、 AP 和 AB的斜率，而在理论上它们对于该点一阶导数的逼近度则分别可从以下泰勒公式的展开式中得知，即由,可见， 对应于式（4-2）和式（4-3）， 它们都截断于 hf(x0)项， 而把 h2项和更高幂次的项全部略去。 换句话说， 就式（4-2）、 式（4-3）而言， 略去余数项所引入的误差将大致和 h 的一次方

6、成正比。,而对于式（4-4）的一阶中心差商表达式则相当于把相应的泰勒公式,截断于 2hf(x0)项， 略去了 h3项以及更高幂次的项。很明显，三种差商表达式中以式（4-4）所示的中心差商的截断误差最小，其误差大致和 h 的二次方成正比。,二阶导数同样可近似为差商的差商，即,这相当于把泰勒公式,截断于 h2f(x)项， 略去了 h4项以及更高幂次的项，其误差亦大致和 h 的二次方成正比。,由此， 仿照式（4-2）和式（4-5）， 偏导数也可近似地用相应的差商来表达。 若设定函数 u (x, y, z)， 当其独立变量 x 得到一个很小的增量x = h 时， 则 x 方向的一阶偏导数可以近似表达为

7、,同样，相应的二阶偏导数可以近似表达为,4.3 差分格式的构造,现以二维静态电、 磁场泊松方程的第一类边值问题为例， 来具体阐明有限差分法的应用。设具有平行平面场特征的电磁场场域 D， 如图 4-2 所示，为一由闭合边界 L 所界定的平面域，其定解问题可表述为,4.3.1 偏微分方程的离散化五点差分格式,通常采用完全有规律的分布方式， 这样在每个离散点上就能得出相同形式的差分方程， 有效地提高解题速度，因而经常采用正方形网 格的剖分方式。现即以这种正方形网格剖分场域 D， 也就是说，用分别与 x、 y 两坐标轴平行的两簇等距（ 步距为 h）网格线来生成正方形网 格， 网格线的交点称为 节点，这

8、样，场域 D 就被离散化为由网格节点构成的离散点的集合。,对于场域内典型的内节点 o (xi，yj)， 如图 4-2 所示， 它与周围相邻的节点 1、 2、 3 和 4构成一个所谓对称的星形。今采用双下标(i，j)的识别方法，设在这些离散节点上的待求位函数 u 的 近 似 值 分 别 记 作 uo = u（i，j）、 u1 = u(i+1，j)、 u2 = u(i，j+1)、 u3 = u(i-1，j) 和u4 = u(i，j-1)， 则参照式(4-7)，二维泊松方程（4-8）可近似离散化表示为,即,此式称为对应于泊松方程的差分方程。如果位函数 u 满足的是拉普拉斯方程（即令式（4-8）中的右

9、端项 F = 0），则差分离散化后所得差分方程是,出现待求函数 u 在点 o(xi，yj)与其四个邻点上的值，故通常称为五点差分格式。,边界条件，对具体问题中可能存在的衔接条件，进行差分离散化处理。,4.3.2 定解条件的离散化各类差分计算格式,对于场域边界上给定的三类边界条件（见 1.7 节）， 由于第二类边界条件可以看作为第三类边界条件的特殊情况，因此，这里只需讨论第一、第三类边界条件的差分离散化处理。,（1） 第一类边界条件的差分离散化,若如图 4-2 点 M 所示， 划分网格时相应的网格节点恰好落在边界 L 上，则只要直接把位函数 u| M L = f(rM)的值赋给该对应的边界节点

10、M 即可。,若划分网格时引入的节点不落在边界 L 上， 则如图 4-3所示， 对于邻近边界的典型节点 o， 由于 h1 h2 h， 这样， o点及其周围相邻的 1、 2、 3 和 4 点构成一个不对称的星形。此时， 可仿照 4.2 节， 采用泰勒公式进行差分离散化处理，即能相当精确地导出关于 o 点的差分计算格式。,应用二元函数的泰勒公式，节点 1 的位函数值 u1 可通过 u0 表示为,同理,以 h 和 h1 分别与以上两式相乘，且相加，然后截断于 h 的二次项，便得关于 的差分表达式为,同理可得,令 h1 =h， h2 =h，代入以上两式，最终再代入给定的泊松方程，即得这类边界情况所对应的

11、差分计算格式为,（2） 第三类边界条件的差分离散化,对此，同样需分两种情况讨论。第一种情况是在边界处引入的相应节点恰好落在边界 L上。 这时，取决于边界 L 在该边界节点处的外法线方向是否与网格线相重合， 对应有不同的差分离散化结果。,当边界 L 在边界节点 o 处的外法向 n 与网格线相重合时，如图 4-4 所示，则问题在于如何用差商近似替代法向导数 。 显然， 最简洁的处理方法是依据式（4-3）， 这样， 第三类边界条件在此情况下的差分计算格式为,当边界 L 在边界节点 o 处的外法向 n 与网格线不重合时，如图 4-5 所示， 显然有,于是， 关于 o 点的差分计算格式是,第二种情况是在

12、边界处引入的相应节点不落在边界 L 上， 这时如图 4-6 所示，可在邻近边界的节点 o 上仍按上述方法列出差分计算格式，只是需引入与节点 o 相关的边界节点 o，取点 o处的外法向 n 作为点 o 处的“外法向 n”， 且近似地认为边界条件中给定的函数f1(ro)和 f2(ro)均在点 o上取值。这样，将式（4-14）中的 f1(ro)和 f2(ro)改记为 f1(ro)和f2(ro)，即得此种情况下关于 o 点的差分计算格式。,应当指出，从实际电、磁场问题的分析需要出发， 如图 4-7 所示， 以通量线（如 E 线）为边界的第二类齐次边界条件是常见的一种情况。这时，边界条件的差分离散化可沿

13、着场域边界外侧安置一排虚设的网格节点， 显然， 对于边界节点 o， 由于该处 ， 故必有 u1 = u3，因此相应于第二类齐次边界条件 的差分计算格式为,4.3.3 不同媒质分界面上边界条件的差分计算格式,当给定的边值问题含有多种媒质时，取决于不同媒质的电磁特性和不同媒质分界面的几何形状， 将对应有类型繁多的差分计算格式，这里仅选取两种典型情况进行分析。,（1） 分界面与网格线相重合的情况,以二维电场问题为例，设分界面 L 与网格线相互重合，如图 4-8 所示。且设在媒质a 中位函数 ua 满足泊松方程，而在媒质b 中位函数 ub 满足拉普拉斯方程。,现若将媒质b 换以媒质a， 则对于 o 点

14、， 据式（4-10）可得,同理，若将媒质a换以媒质b， 则对于 o 点， 据式（4-11）可得,但实际上 ua1和 ub3是虚设的电位，所以应利用分界面上场量遵循的边界条件式（1-66）和式（1-69），把它们从以上两式中消去。,首先， 由式（1-66）得出分界面上电位的连续性，即,其次，假设在分界面上自由电荷的面密度 = 0， 则由式（1-69）有,以差分格式表示，即为,将a 乘以式（4-16）与b 乘以式（4-17）后相加， 代入由式（4-18）和式（4-19）所给定的边界条件， 并令 K =a/b，便得待求的两种不同媒质分界面上边界条件的差分计算格式为,（2） 分界面对于网格呈对角线形态

15、的情况,此时，差分计算格式的推导及其处理方法与上类同， 但为提高差分离散化的逼近度， 尚需引入 M、 N 两个辅助节点（ M、 N二点分别是线段14和23的中点）， 如图 4-9 所示。对于节点 o， 如同前述，当媒质依次代换时， 相应的五点差分格式分别与式（4-16）和式（4-17）相同。依据分界面上的边界条件，现应有,注意到在以上各式中 ua1、 ua4、 uaM、 ub2、 ub3和 ubN都是虚设电位值，但应用线性插值，它们可由以下方程相互关联：,因此，由实际存在的电位值,可以消去所有虚设电位值，得出关于这类边界条件的差分计算格式,4.3.4 对称线的差分计算格式,在实际分析电、 磁场

16、分布时， 经常可观察到场分布的对称性， 因此，在数值计算中计及场的对称线条件，即可缩小分析计算的场域，从而在对计算机存贮容量要求不变的情况下，可获得更为理想的数值解。,设如图 4-10 所示， AA线为二维泊松场的对称线。此时，对位于对称线上的任一节点 o， 由 式（4-10）， 并依据场的对称性，即有 u1 = u3，因此相应的差分计算格式为,4.4 差分方程组的求解,综上所述，对场域 D 内各个节点（包括所有场域内节点和有关的边界节点）逐一列出对应的差分计算格式，即构成以这些离散节点上的位函数 u 为待求量的差分方程组（代数方程组）。仔细分析所得的差分方程组，不难看出，该方程组的系数一般都

17、是有规律的， 且各个方程都很简单，包含的项数不多（取决于前述对称或不对称的所谓星形离散结构， 每个方程待求量的项数最多不超过 5 项）。 因此， 在第二章所述的众多代数解法中，对于有限差分法，通常都采用迭代法，这是因为用计算程序来实现迭代时， 需要用到哪些系数就算出哪些系数，不需用时不保留，这样可显著降低对计算机存贮容量的需求。,在迭代法的应用中，为加速迭代解的收敛速度，通常采用的是逐次超松弛迭代法。按图 4-11 所示的对称星形离散模式，对应于泊松差分方程（4-10）， 若采用早期的高斯赛德尔迭代法（规定迭代运算顺序是： 从左下角开始做起， 即 i 小的先做； 对固定的 i，j 小的先做。）

18、， 则关于节点 o 迭代到第（n +1）次时的近似值，应由如下迭代公式算得,而为加速迭代解的收敛， 构成超松弛迭代公式的原则是： 并不将由上式所算得的结果作为u(i，j)的第（n +1）次近似值，而仅把它视为一中间结果 然后作加权平均处理，即令,式中， 称为加速收敛的松弛因子。 很明显， 上式就是 2.5 节中已经给出的一般计算公式（2-14）对应于本方法的具体表达式。正如前已指出的， 超松弛迭代法的 取值范围是 1 2， 当 =1 时， 式（4-28）即归结为高斯赛德尔迭代法的迭代公式（4-27）； 当 2 时，迭代过程将不收敛而发散。最佳收敛因子的取值随问题和离散化的情况而异。对于第一类边

19、值问题，若一正方形场域由正方形网格剖分（每边节点数为 p +1）， 则最佳收敛因子opt可按下式计算,若一矩形场域由边长为 h 的正方形网格剖分（ 设两边分别为 ph 和 qh， 且 p、q 通常要大于15），则相应的最佳收敛因子为,在更一般的情况下， opt只能凭经验取值。 值得指出， 在 2.5 节中，介绍了加速收敛的松弛因子 作自适应估计的方法， 这为解决一般性的需要提供了优化加速收敛因子选择的数学工具，然而，这时不仅首先必须形成差分方程组所对应的系数阵，而且相继需要构造系数矩阵元素的存贮技术（如 2.5.3 节所阐述的非零元素存贮技术）。 换句话说，应用数学上的高要求导致了求解过程的复

20、杂化。,应当注意，在迭代运算前，恰当地给定各内点的初值（即所谓零次近似值），也是加速收敛速度的一个有效途径。,在超松弛迭代法的应用中，还必须涉及迭代解收敛程度的检验问题。对此，通常的处理方法是：以所有内点上相邻两次迭代解的绝对误差或相对误差不大于指定的误差范围，作为检查迭代解收敛程度的依据。,4.5 场强与电、磁积分量的计算,通过上述差分方程组的求解，在获得待求位函数 u(x,y)的数值解后， 往往还需求场中的场强分布，以及其他有关的积分特性（如磁通量和磁导、电导、电容等磁路及电路参数等）。现以二维平行平面场为例导出关于这些物理量和参数的差分计算公式，推导中设场域由正方形网格予以剖分。,4.5

21、.1 场强的差分计算公式,基于 1.6 节的阐述，在静态二维场中， 电场强度 E、 磁场强度 H 或磁感应强度 B 和它们对应的位函数之间的关系可用差商分别表示为,式中， M、 Mm和 MA 分别为电位、标量磁位和向量磁位函数的标度，定义为相应位函数的实际值与相对值之比。例如在 4.6 节例 4-1 中， 采用了1 =10 的相对电位值， 若1 的实际值为 150V， 则计算电场强度时引入的电位函数标度应该是 M=150V/10 =15V。 但若计算时，位函数直接采用实际值，则 M=1。,如需计算边界上的场强，由于按式（4-31） 式（4-33）中所取的位函数值通常是在相距为h 而非 2h 的

22、两点上的值，因此所得结果实际上并不是边界处的场强， 而应该是与边界相邻的网格边和边界的中间点上的场强值。例如，图 4-15d 中边界点 S 上的电场强度即可表示为,显然，只有当网格的步距 h 足够小时， 上式计算结果才有可能逼近边界点 S 上实际的场强值。,4.5.2 通量与参数的差分计算式,无论是静电场、恒定电流场或恒定磁场，其通量 可一般性地表示为,式中， K 是相应媒质的宏观特征参数（、 或）， 而 a 则为上述各类电、 磁场的相关场量（即相应场强 E、E或 H）。在求得场中各点场强的基础上，这一通量积分值可以近似地表示成,式中， n 表示被积面积被网格剖分所得小块面积的总数； Si 表

23、示其中某一小块面积； aav(i)表示在 Si 中所取的场强 ai 的平均值， 并且 ai 的方向应与小面积 Si 的法线方向相一致。 这样， 在通量的差分计算式（4-35）的基础上，所分析的静电场中的电容 C、 恒定电流场中的电导 G 或恒定磁场中的磁导 等电路或磁路参数 P 就可按下式计算：,式中， U 表示限定分析区域的边界面间的电位差或磁位差。,例 4-2 二维平行平面电流场的计算。,在导电纸模拟的实验研究中，制备了如图 4-17 所示的两维电流场模型，其中两种导电媒质的电导率分别为1 和2， 它们在场域的对角线 L上接合。 电极间外施电压 10V。试求该电流场模型中两维电流场分布。,

24、4.7 等值点的寻求与描绘,在电磁场分布的研究中，为了形象化的分析需要，通常需要通过数值计算的后处理，描绘出场分布的可视化图形，从而可定性乃至定量地讨论场分布的规律性。常见的场分布图形为电场中的等位面（线）、 磁场中的等磁位面（线）以及磁感应强度 B 线的分布等。 应再次指出，诚如 3.5 节的讨论， 在具有平行平面场或轴对称场特征的前提下， 借助于向量磁位 A(Az =const. 或A =const.)即可方便地描绘出相应磁场的 B 线分布。这些由相应的位函数数值相等的点所形成的曲面（线），称为等值面（线），其一般方程为,对少量电磁场问题，上式可由解析表达式给出， 利用该表达式就可以直接绘

25、制场分布图形。对大量的工程电磁场问题来说，则必须有赖于电磁场的数值解，通过插值法来寻求对应于给定位值的等值点。,4.7.1 等值点的寻求,以平行平面电场中等值点的寻求为例， 当由有限差分法算出各网格节点的电位值后， 可以利用线性插值关系来求得指定电位值的等值点坐标 (x, y)。具体方法和步骤如下：,（1） 给出等值线的指定电位值 Veq；,（2） 判断相应的网格线是否与位值等于给定的 Veq的等值线相交。 如图 4-21 所示， 设某个正方形网格的四顶点坐标分别为 A(i, j)、 B(i, j+1)、 C(i+1, j+1)和 D(i+1，j)， 现首先判断网格线 AB 是否与指定位值的等

26、值线相交。 显然， 若点 A 和点 B 的电位值 V(i，j)和 V(i，j+1)与指定位值之间满足下列不等式：,则指定的等值线必与网格线 AB 相交，换句话说，在 AB 线上有相应的交点存在。,(3） 按线性插值关系，确定上述所得等值点的坐标。在网格线分别沿 x， y 坐标轴取向的前提条件下，按线性插值公式,即可求得上述交点（有指定位值的等值点）的坐标为,（4） 同理，继续搜索沿 x 方向网格线 AD 上是否存在待求的等值点。一旦存在，则其计算关系式可类同推得为,至此，对各个网孔分别在相应的 x 和 y 方向的两网格线上搜索对应的等值点；为使所得等值点形成有序的排列，这里，还运用所谓冒泡法对

27、选定的某个坐标方向实现等值点的排序处理。冒泡法的思路在于将相邻的两个数值进行比较，将数值小的一个调迁到前一位置。,4.7.2 等值面（线）的绘制,对应于场分布（等值面或等值线）图形描绘的需求，在寻找出各组等值点的分布后，由所构成的数据文件，即可借助于各类绘图软件，例如 Math、 Tech* Graph* Pad 等， 以及如常用软件 MATLAB 等，完成等值面（线）的绘制。书在 5.4.5 节对等值线的绘制，基于三角元剖分，进行了系统的展述。,4.8 时域有限差分法,近代技术的发展，使复杂的高频电磁系统的分析与综合，以及高频电磁场与复杂目标相互作用的分析和计算，成为重要的研究课题。这些研究

28、课题以高频电磁场的传输、辐射、散射、和透入问题为主线，反映了现代通信、 雷达、 物探、 电磁防护、 电磁兼容、 医疗诊断、战略防御以及工农业生产和日常生活等领域多方面的需求。正是在众多分析任务与目标的推动下，时域有限差分法历经 20 余年的发展，以其直接的时域计算模式、广泛的适用性、较经济的存贮空间和计算时间、程序的通用性与简明、直观等特点，从传统的有限差分法中脱颖而出，成为在上述一系列研究课题中广泛应用的数值计算方法。本节即在于概述时域有限差分法（FDTD）的基本应用原理。,1966 年 KaneS.Yee 提出了后被称为 Yee 氏网格的空间离散方式（见图 4-22）。 这一合理的网格体系

29、的特点是，电场和磁场各分量在空间的取值点被交叉地放置，从而在直角坐标系下每个坐标平面上相应的电场分量的四周由磁场分量环绕，而相应的磁场分量的四周则由电场分量环绕。这样的网格空间配置符合法拉第电磁感应定律和安培环路定律的要求。例如对应于图 4-22b 中环绕点 o(xo, yo,zo)的环量,4.8.1 Yee 氏网格,式中，在相应的元路径x 或y 上求积时， 对应场量 Ex 或 Ey 被看作为常量， 且分别等于元路径中点处的 Ex 或 Ey 值。从而通过应用二元函数的泰勒公式， 并截断于一阶偏导数项，可得,将以上关于 Ex1、 Ey2、 Ex3和 Ey4的近似表达式代入式（4-39）， 即有,

30、而依据法拉第电磁感应定律应有,综合式（4-40）和式（4-41），显然满足麦克斯韦方程组中的旋度方程 , 即,由此可见， Yee 氏网格体系反映了实际物理模型中电场和磁场互为因果的物理本质，即满足麦克斯韦方程组的两个旋度方程，因而也就符合电磁波在空间传播的规律性。此外，它也满足不同介质分界面上场的切向分量连续的物理条件。 显然， Yee 氏网格为在四维空间中合理地离散六个未知场量，建立具有高精度的差分计算格式，奠定了理想的离散化空间的应用基础。,4.8.2 旋度方程的差分格式,当场域由 Yee 氏网格离散后，空间步长分别为 x、 y 和 z； 时间步长记为 t， 以 n表示时间步长的“个数”，

31、并标记于右上角，因而场分量 F（x，y，z，t）的四维空间离散表示法为,在 1.4.1 节中，已经指出，麦克斯韦方程组中两旋度方程是基本的，这是电磁场问题研究的出发点。,应指出，为保证 FDTD 计算稳定性，时间离散的步长与空间离散步长间应满足一定的关系。经分析表明，时间步长可选为电磁波传播一个空间步长所需时间的一半。,现应用中心差商近似替代该场分量对空间、时间的偏导数，即,在无源、均匀且各向同性的线性介质中，麦克斯韦方程组的两旋度方程分别为,以上两式在直角坐标系下的展开式，分别为,应用式（4-43）和式（4-4）的差商近似关系式， 式（4-47a）对应的差分计算格式为,同 理， 在 (n +

32、1/2)时 间 步， 对 (i， j+1/2，k )点 的 Ey； 在 (n +1/2)时 间 步， 对(i， j， k +1/2)点的 Ez， 可得与式（4-47b）和式（4-47c）分别对应的完全类似的差分格式。,对于相应的第二旋度方程中的磁场分量的差分格式，由方程的对称性，可类比求得。应注意的是， 因在 Ex、 Ey 和 Ez 差分格式中磁场值取于 (n +1/2)时间步， 故下式中磁场取值均应取自 (n +1/2)时间步或 (n -1/2)时间步，以保证取值的时间步差为一个整时间步， 从而保证下式中电场分量取值时间与前面的电场分量取值时间相同。这样，将为未知量的存贮和计算带来很多方便。

33、因此，对于式（4-48a）应有,同理， 在 n 时间步， 对( i+1/2， j， k +1/2)点的 Hy； 在 n 时间步， 对 (i+1/2， j+1/2，k )点的 Hz， 可得与式（4-48b）和式（4-48c）分别对应的完全类似的差分格式。,可以看出，任一网格点上的电场值只与它上一时间步的电场值及四周环绕它的磁场值相关；同样，任一网格点上的磁场值也只与它上一时间步的磁场值及四周环绕它的电场值相关。此外，媒质参数、均为空间坐标的函数， 故 FDTD 易于处理非均匀和各向异性媒质的问题。,4.8.3 解的数值稳定性,在时域有限差分法中，时间增量 t 和空间增量 x、 y 和 z 不是相

34、互独立的， 它们之间必须满足一定的关系， 否则， 将出现算法上的不稳定性。可以证明10 ， 在三维情况下， FDTD 数值稳定条件为,式中， V 是电磁波的传播速度。,如果x =y =z =l 时， 则 ， 一般取,若在三个坐标轴方向上的空间步长是可变的， 那么， 应该先取每一坐标方向上的最小步长，然后三者中再选最小值，即,4.8.4 边界条件,关于边界条件包含三种类型：有界场域的边界条件；不同介质分界面上的边界条件和无界场域截断边界上的吸收边界条件。,对于切向电场或法向磁场为零的边界，如 5.7.2 节所述的波导场问题， 此时在波导壁边界上， 只要使切向电场或法向磁场所在的网格节点落在边界上

35、，则对于 TM 波， 应令边界节点上的 Et = 0； 对于 TE 波，由于磁场的法向导数为零， 可虚拟与边界面相距半个步长的网格， 并令节点上的值对于边界面呈偶对称分布，其分析方法与 4.3.2 节中所表述的边界条件处理方法类似。,（1） 有界场域的边界条件,（2） 不同介质分界面上的边界条件,设典型的两种理想介质分界面如图 4-23 所示， 其分界面 S 与坐标面 XOZ 相平行。在两种介质区域中， 式（4-47a）应分别为,在介质分界面 S 上，磁场分量连续，即 Hz1 Hz2 = Hz， Hy1 Hy2 = Hy（理想介质的磁导率为0）； 电场的切向分量连续， 即 Ex1 Ex2 =

36、Ex。 将上述关系代入式（4-54）和式（4-55）， 然后两式相加，得到,式中， =1/2(1 +2)。可以看出， 式（4-56）与式（4-47a）的表达形式完全相同。因此，在建立所述理想介质分界面上边界条件所对应的差分格式时，只要将式（4-49）中的 用代替即可。同理，可相仿处理另两个典型方位的不同介质分界面上的边界条件。,（3） 吸收边界条件,在无限长波导、电磁辐射、电磁散射等问题中，场域延展至无限远处，电磁波没有反射效应。可以设想，如果对所感兴趣的区域划定一个边界，使电磁波在该边界上没有明显地反射，则就不会因该界面的引入导致内部空间场分布发生畸变。换句话说，在该限定的有限空间中的计算结

37、果将与原无限空间中的计算结果等效。具有这种功能的边界条件称为吸收边界条件或辐射边界条件。这里简单介绍 Mur 提出的如下吸收边界条件。,应指出，在时变电磁场中，如果在介质分界面上场量满足切向分量连续的条件，则该场量将自动满足法向分量连续的条件。也就是说，在讨论吸收边界条件时，只需考虑截断面上场量的切向分量即可。,以 F 表征场分量，可以证明，对于三维问题，在截断面 x = 0 和 x = xm 上的吸收边界条件分别为,式中， 为在理想介质中的电磁波波速。 其他两个方向上的吸收边界条件可以类同推得。,对于二维问题，吸收边界条件分别为,另一方向上的吸收边界条件可以类同推得。,对于一维问题，吸收边界

38、条件分别为,不难看出，吸收边界条件不能直接应用前述的差分格式式（4-43）， 这是因为中心差商格式是由所论网格节点四周相邻半个网格尺寸的场值决定的，而在截断边界处，边界以外的半个网格节点上的值是未知的。,现从一维吸收边界条件入手阐述。设均匀平面波在无界空间沿 x 轴方向传播，研究区域为x0，xm区间。 为此， 在 x = x0 和 x = xm 处分别设置截断点， 如图 4-24 所示。 在边界 x = x0 处， 引入辅助场量 F(1/2)， 它位于节点 0 与1 之间， 从而式（4-61）的中心差商格式为,对 Fn +1/2(1) 取时间相邻的节点平均值， 即 Fn +1/2 =1/2( Fn +1+ Fn）对 F(1/2)取空间相邻节点的平均值，即 F(1/2)=1/2 F(1)+ F(0) 。由此代入式（4-63）， 便得 x = x0 处一维吸收边界条件的差分格式为,同理， 对于式（4-62）， 在边界 x = xm 处的一维吸收边界条件的差分格式为,对于二维情况，由式（4-59）可得,对于 x = xm 边界面，以及另一坐标轴方向上的边界条件， 乃至三维情况下的边界条件， 均可由相似的推导过程得出。,