《计算电磁学》第五讲

2020年4月11日9时50分第五讲吸收边界条件Dr.PingDu(杜平)E-mail:[e-mail]ology2020年4月11日9时50分为何要用吸收边界条件？由于计算机内存空间有限，让我们的分析区域不能无限大，必须在某处截断。另外，即使内存空间很大，也必须在某处截断。因为大的分析空间就意味着巨大的计算时间。从节省计算时间的角度考虑，也必须截断。2020年4月11日9时50分有哪些吸收边界条件？©基于Sommerfeld辐射条件的Bayliss-Turkel吸收边界条件；©基于单向波动方程的Engquist-Majda吸收边界条件(1977年)；©Mur吸收边界条件(1981)；©Trefethen-Halpern近似展开(1986)©Higdon算子；2020年4月11日9时50分©利用差值技术的廖氏吸收边界条件廖氏吸收边界条件比Mur二阶吸收边界条件在网格外边界引起的反射要小一个数量级(20dB)，对外向波传播角度或数值色散均不敏感，并且在矩形计算区域的角点处也易于实现。M.Moghaddam,R.L.Wagner和W.C.Chew(周永祖)曾指出，如果采用单精度计算，可能导致使用廖氏吸收边界条件的FDTD算法不稳定，而采用双精度则可改善稳定性。©Mei-Fang超吸收边界条件(1992)；©PML完全匹配层2020年4月11日9时50分Engquist-Majda吸收边界条件考虑二维情形时的齐次波动方程，222222210UUUxyct(5-1)其中U为标量场分量，c为波的相速度。定义偏微分算子2222222222211xytLDDDxyctc(5-2)于是方程(5-1)可写为0LU(5-3)2020年4月11日9时50分算子L还可以通过因式分解写为0LLU(5-4)其中，21txDLDsc21txDLDsc/ytDsDc(5-5)(5-6)(5-7)2020年4月11日9时50分在网格边界，如x=0处，将算子作用于波函数将完全吸收以任意角度入射到边界的平面波，即将L0LU(5-8)用于图1中的边界x=0，可构成一个准确的解析吸收边界条件。它将吸收来自区域Ω内的波。iIncidentwaveReflectionwaveaboxy图1二维吸收边界条件Fig.1.Two-dimensionalABC.2020年4月11日9时50分相似地，算子作用于波函数，将构成x=a处的准确吸收边界条件。L对式(5-5)和(5-6)中根式的处理，可以用Taylor级数展开。将21s在s=0附近展开为Taylor级数，2211...2ss(5-9)当s很小时，只取一项，211s(5-10)2020年4月11日9时50分将(5-10)代入(5-5)中，有txDLDc(5-11)将其代入(5-8)，可得10UUxct这就是所分析区域左侧边界x=0的一阶近似吸收边界条件。(5-12)2020年4月11日9时50分将(5-9)中的级数取两项，有22112ss(5-13)Substitution(5-13)in(5-5)yields2221111222yttxxtytxtcDDDLDsDccDcDDDcD(5-14)2020年4月11日9时50分Substituting(5-14)in(5-8)has22222102UUcUxtcty这就是所分析区域左侧边界x=0的二阶近似吸收边界条件。对于图1中的其他边界，相应的二阶近似解析吸收边界条件为22222102UUcUxtcty，x=a边界（右边界）(5-15)(5-16)2020年4月11日9时50分22222102UUcUytctx，y=0边界（下边界）22222102UUcUytctx，y=b边界（上边界）(5-17)(5-18)考虑三维情形时的齐次波动方程，22222222210UUUUxyzct(5-19)2020年4月11日9时50分此时，偏微分算子2222222222222211xyztLDDDDxyzctcLL将L分解为和,得到与(5)和(6)相同的准确吸收边界条件算子。不同的是，s为(5-20)22//yzttDDsDcDc(5-21)2020年4月11日9时50分L算子作用于波函数U，将在网格左边界x=0处准确地吸收以任意角度入射到边界的平面波。利用Taylor级数近似展开式(5-10)，可得到x=0处的一阶吸收边界条件，其形式与(5-12)相同。利用Taylor级数近似展开式(5-13)，可得到x=0处的二阶吸收边界条件。其表达式为22022ytzxttcDDcDDUcDD(5-22)2020年4月11日9时50分两边同乘以，得tD22222221022UUcUcUxtctyz(5-23)当s很小时，(5-23)是准确吸收边界条件的很好近似。LU对于其他网格边界相应的二阶近似解析吸收边界条件为22222221022UUcUcUxtctyz，x=a边界(5-24)2020年4月11日9时50分，y=0边界，y=b边界，z=0边界22222221022UUcUcUytctxz22222221022UUcUcUytctxz(5-25)(5-26)22222221022UUcUcUztctxy(5-27)22222221022UUcUcUztctxy，z=h边界(5-28)对于矢量Maxwell方程的FDTD仿真，近似吸收边界条件(5-23)-(5-28)中的U表示位于网格边界上的E和H的各个切向分量。2020年4月11日9时50分Mur差分格式对于上述一阶、二阶近似解析吸收边界条件，Mur提出了一种简单有效的差分数值算法。利用它们来截断FDTD仿真区域，总体虚假反射在1%-5%。以一阶情形，x=0边界为例。在12xx(1/2)tnt处、时刻用中点差分来代替式(5-12)中的偏微分，得112211,()122(1)(0)nnUxntUUxx(5-29)2020年4月11日9时50分111,()111122()()22nnUxntUUctct(5-30)其中半网格点和半时间步长时刻的值，可用下列二阶精度的平均公式计算1121()()()2nnnUmUmUm11()(1)()22nnnUmUmUm(5-31)(5-32)2020年4月11日9时50分Substituting(5-31),(5-32)in(5-29),(5-30)and(5-12),afterthemanipulation,gets11(0)(1)(1)(0)nnnnctxUUUUctx(5-33)这就是Mur一阶吸收边界条件(Mur’s1storderABC)。下面推导Mur二阶吸收边界条件。仍旧以x=0边界为例，如图2所示。2020年4月11日9时50分(0,1)Wj(0,)Wj(0,1)Wj1(,)2WjX(1,1)Wj(1,)Wj(1,1)Wj(1)yjyyjy(1)yjy0xxx图2Mur吸收边界条件Fig.2.Mur’s2ndABC.图中(0,)nWj表示位于x=0网格边界的E或H的切向分量。2020年4月11日9时50分具体为，在距离网格边界半步长的辅助网格点1(,)2j，对1(,)2nWj分量,将(5-23)中的偏微分用中心差分来代替。将式中关于x、t的偏导用中心差分格式写为112111(,)(,)(,)2221111(1,)(0,)(1,)(0,)1212nnnjjjnnnnjjjj(5-34)2020年4月11日9时50分将该式对时间t的偏导数写为相邻两点(0,j)和(1,j)处对时间偏导的平均，其表达式为2222221(,)(0,)(1,)21111(0,)(0,)(0,)(1,)(1,)(1,)22122212()()nnnjjjnnnnnnjjjjjj(5-35)2020年4月11日9时50分(3)将该式对y的偏导数写为相邻两点(0,j)和(1,j)处对y偏导的平均，其表达式为2222221(,)(0,)(1,)2(0,1)(0,)(0,1)(1,1)(1,)(1,1)22122212()()nnnjjjnnnnnnjjjjjj(5-36)(4)Substituting(5-34)-(5-36)in(5-23),onecanobtain1(0,)njW.Thus,the2ndABCatx=0boundarycanbegot,whichis2020年4月11日9时50分1111(0,,)(1,)(1,)(0,)(1,)(0,)22(0,1)(0,)(0,1)(1,1)(1,)(1,1)2()22()()2nnnnjkjjjnnjjnnnjjjnnnjjjctx(5-37)2020年4月11日9时50分对于三维情形，考虑x=0边界。从式(5-23)出发，此时前面的网格图2位于zkz网格平面，在距离边界半个步长的辅助网格点处，对分量1(,,)2jk1(,,)2nWjk用中点差分代替(5-23)中的偏导运算。2Wxt22Wt22Wyzkz偏导、和与式(5-34)-(5-36)相同。只是在平面进行计算。2020年4月11日9时50分偏导可表达为在(0,j,k)和(1,j,k)处对z的偏导的平均值，其为22Wz2222221(,,)(0,,)(1,,)2(0,,1)(0,,)(0,,1)(1,,1)(1,,)(1,,1)22122212()()nnnjkjkjknnnnnnjkjkjkjkjkjk(5-38)将这些差分表达式代入式(5-23)，解出就得到三维情形下W分量在x=0网格边界出的二阶吸收边界条件10,,njkW2020年4月11日9时50分1111(0,,)(1,,)(1,,)(0,,)(1,,)(0,,)22(0,1,)(0,,)(0,1,)(1,1,)(1,,)(1,1,)222()22()()2()2()nnnnjkjkjkjknnjkjknnnjkjkjknnnjkjkjkctx(0,,1)(0,,)(0,,1)(1,,1)(1,,)(1,,1)2()2nnnjkjkjknnnjkjkjk(5-39)2020年4月11日9时50分问题：边界反射系数与入射角的关系。葛德彪、闫玉波《电磁波时域有限差分方法》一书给出了Mur近似边界条件的反射系数