《计算电磁学》第五讲

1、2022年4月23日0时52分第五讲 吸收边界条件Dr. Ping Du(杜平)E-mail: HeFei University of Technology2022年4月23日0时52分为何要用吸收边界条件？ 由于计算机内存空间有限，让我们的分析区域不能无限大，必须在某处截断。另外，即使内存空间很大，也必须在某处截断。因为大的分析空间就意味着巨大的计算时间。从节省计算时间的角度考虑，也必须截断。 2022年4月23日0时52分有哪些吸收边界条件？ 基于Sommerfeld辐射条件的Bayliss-Turkel吸收边界条件； 基于单向波动方程的Engquist-Majda吸收边界条件(1977年

2、)； Mur吸收边界条件(1981) ； Trefethen-Halpern近似展开(1986) Higdon算子；2022年4月23日0时52分 利用差值技术的廖氏吸收边界条件 廖氏吸收边界条件比Mur二阶吸收边界条件在网格外边界引起的反射要小一个数量级(20dB)，对外向波传播角度或数值色散均不敏感，并且在矩形计算区域的角点处也易于实现。 M. Moghaddam, R. L. Wagner和W. C. Chew (周永祖) 曾指出，如果采用单精度计算，可能导致使用廖氏吸收边界条件的FDTD算法不稳定，而采用双精度则可改善稳定性。 Mei-Fang超吸收边界条件(1992) ； PML完全

3、匹配层2022年4月23日0时52分n Engquist-Majda吸收边界条件 考虑二维情形时的齐次波动方程， 222222210UUUxyct(5-1)其中U为标量场分量，c为波的相速度。定义偏微分算子 2222222222211xytLDDDxyctc(5-2)于是方程(5-1)可写为 0LU (5-3)2022年4月23日0时52分算子L还可以通过因式分解写为 0L L U(5-4)其中， 21txDLDsc21txDLDsc/ytDsDc(5-5)(5-6)(5-7)2022年4月23日0时52分在网格边界，如x=0处，将算子 作用于波函数将完全吸收以任意角度 入射到边界的平面波，即

4、将 L0L U(5-8)用于图1中的边界x=0，可构成一个准确的解析吸收边界条件。它将吸收来自区域内的波。 图1 二维吸收边界条件Fig.1. Two-dimensional ABC. 2022年4月23日0时52分相似地，算子 作用于波函数，将构成x=a处的准确吸收边界条件。 L对式(5-5)和(5-6)中根式的处理，可以用Taylor级数展开。 将21 s在s=0附近展开为Taylor级数， 2211.2ss (5-9)当s很小时，只取一项， 211s(5-10)2022年4月23日0时52分将(5-10)代入(5-5)中，有 txDLDc(5-11)将其代入(5-8)，可得10UUxct

5、这就是所分析区域左侧边界x=0的一阶近似吸收边界条件。 (5-12)2022年4月23日0时52分将(5-9)中的级数取两项，有 22112ss (5-13)Substitution (5-13) in (5-5) yields 2221111222yttxxtytxtcDDDLDsDccDcDDDcD(5-14)2022年4月23日0时52分Substituting (5-14) in (5-8) has 22222102UUcUx tcty 这就是所分析区域左侧边界x=0的二阶近似吸收边界条件。 对于图1中的其他边界，相应的二阶近似解析吸收边界条件为 22222102UUcUx tcty

6、，x=a边界（右边界） (5-15)(5-16)2022年4月23日0时52分22222102UUcUy tctx ，y=0边界（下边界） 22222102UUcUy tctx ，y=b边界（上边界） (5-17)(5-18)考虑三维情形时的齐次波动方程， 22222222210UUUUxyzct(5-19)2022年4月23日0时52分此时，偏微分算子 2222222222222211xyztLDDDDxyzctcLL将L分解为和,得到与(5)和(6)相同的准确吸收边界条件算子。不同的是，s为(5-20)22/yzttDDsDcDc(5-21)2022年4月23日0时52分L算子作用于波函数

7、U，将在网格左边界x=0处准确地吸收以任意角度入射到边界的平面波。 利用Taylor级数近似展开式(5-10)，可得到x=0处的一阶吸收边界条件，其形式与(5-12)相同。 利用Taylor级数近似展开式(5-13)，可得到x=0处的二阶吸收边界条件。其表达式为 22022ytzxttcDDcDDUcDD(5-22)2022年4月23日0时52分两边同乘以 ，得 tD22222221022UUcUcUx tctyz (5-23)当s很小时，(5-23)是准确吸收边界条件 的很好近似。 L U对于其他网格边界相应的二阶近似解析吸收边界条件为22222221022UUcUcUx tctyz ，x=

8、a边界 (5-24)2022年4月23日0时52分，y=0边界 ，y=b边界 ，z=0边界 22222221022UUcUcUy tctxz 22222221022UUcUcUy tctxz (5-25)(5-26)22222221022UUcUcUz tctxy (5-27)22222221022UUcUcUz tctxy ，z=h边界 (5-28) 对于矢量Maxwell方程的FDTD仿真，近似吸收边界条件(5-23)-(5-28)中的U表示位于网格边界上的E和H的各个切向分量。 2022年4月23日0时52分n Mur差分格式 对于上述一阶、二阶近似解析吸收边界条件，Mur提出了一种简单

9、有效的差分数值算法。利用它们来截断FDTD仿真区域，总体虚假反射在1%-5%。 以一阶情形，x=0边界为例。在 12xx(1/ 2)tnt处、 时刻 用中点差分来代替式(5-12)中的偏微分，得 112211,()122(1)(0)nnUx ntUUxx(5-29)2022年4月23日0时52分111,()111122( )( )22nnUx ntUUctc t(5-30)其中半网格点和半时间步长时刻的值，可用下列二阶精度的平均公式计算 1121( )( )( )2nnnUmUmUm11()(1)( )22nnnUmUmUm(5-31)(5-32)2022年4月23日0时52分Substitu

10、ting (5-31), (5-32) in (5-29), (5-30) and (5-12), after the manipulation, gets 11(0)(1)(1)(0)nnnnc txUUUUc tx (5-33)这就是Mur一阶吸收边界条件(Murs 1st order ABC)。 下面推导Mur二阶吸收边界条件。仍旧以x=0边界为例，如图2所示。 2022年4月23日0时52分图2 Mur吸收边界条件 Fig. 2. Murs 2nd ABC. 图中(0, )nWj表示位于x=0网格边界的E或H的切向分量。 2022年4月23日0时52分具体为，在距离网格边界半步长的辅助

11、网格点 1( , )2j，对1( , )2nWj分量,将(5-23)中的偏微分用中心差分来代替。 将式中关于x、t的偏导用中心差分格式写为112111( , )( , )( , )2221111(1, )(0, )(1, )(0, )1212nnnjjjnnnnjjjjWWWx ttxxWWWWtxx (5-34)2022年4月23日0时52分 将该式对时间t的偏导数写为相邻两点(0, j)和(1, j)处对时间偏导的平均， 其表达式为2222221( , )(0, )(1, )21111(0, )(0, )(0, )(1, )(1, )(1, )22122212()()nnnjjjnnnnn

12、njjjjjjWWWtttWWWWWWtt(5-35)2022年4月23日0时52分(3) 将该式对y的偏导数写为相邻两点(0, j)和(1, j)处对y偏导的平均， 其表达式为 2222221( , )(0, )(1, )2(0,1)(0, )(0,1)(1,1)(1, )(1,1)22122212()()nnnjjjnnnnnnjjjjjjWWWyyyWWWWWWyy(5-36)(4) Substituting (5-34)-(5-36) in (5-23), one can obtain 1(0, )njW.Thus, the 2nd ABC at x=0 boundary can be

13、 got, which is 2022年4月23日0时52分1111(0, , )(1, )(1, )(0, )(1, )(0, )22(0,1)(0, )(0,1)(1,1)(1, )(1,1)2()22() ()2nnnnj kjjjnnjjnnnjjjnnnjjjc txWWWWc txxWWc txc txWWWyc txWWW (5-37)2022年4月23日0时52分对于三维情形，考虑x=0边界。 从式(5-23)出发，此时前面的网格图2位于 zk z 网格平面， 在距离边界半个步长的辅助网格点 处，对分量 1( , , )2j k1( , , )2nWj k用中点差分代替(5-2

14、3)中的偏导运算。 2Wx t 22Wt22Wyzk z 偏导、和与式(5-34)-(5-36)相同。只是在平面进行计算。 2022年4月23日0时52分偏导 可表达为在(0, j, k)和(1, j, k)处对z的偏导的平均值，其为 22Wz2222221( , , )(0, , )(1, , )2(0, ,1)(0, , )(0, ,1)(1, ,1)(1, , )(1, ,1)22122212()()nnnj kj kj knnnnnnj kj kj kj kj kj kWWWzzzWWWWWWzz(5-38) 将这些差分表达式代入式(5-23)，解出 就得到三维情形下W分量在x=0网格

15、边界出的二阶吸收边界条件 10, ,nj kW2022年4月23日0时52分1111(0, , )(1, , )(1, , )(0, , )(1, , )(0, , )22(0,1, )(0, , )(0,1, )(1,1, )(1, , )(1,1, )222()22() ()2()2()nnnnj kj kj kj knnj kj knnnjkj kjknnnjkj kjkc txWWWWc txxWWc txc txWWWyc txWWWc txz (0, ,1)(0, , )(0, ,1)(1, ,1)(1, , )(1, ,1)2()2nnnj kj kj knnnj kj kj k

16、WWWc txWWW (5-39)2022年4月23日0时52分问题：边界反射系数与入射角的关系。 葛德彪、闫玉波电磁波时域有限差分方法一书给出了Mur近似边界条件的反射系数。cos() 1cos() 1iiR(5-40)Mur二阶近似边界条件时的反射系数：Mur一阶近似边界条件时的反射系数：22cos() 1 sin ()/ 2cos() 1 sin ()/ 2iiiiR (5-41)其中 为入射角i2022年4月23日0时52分反射系数的绝对值与入射角的关系，如图3所示。. 010203040506070809001i( o )|R|1st Order Approximation 2nd Order Approximation图3 反射系数的绝对值与入射角的关系Fig. 3. Relation between the reflection coefficient and the incident angle. 2022年4月23日0时52分n 完全匹配层(Perfect Matched Layer, PML)1994年，Berenge