（计算数学专业论文）电磁波fdtd方法的若干研究.pdf

摘要 摘要 本文主要是关于电磁波时域有限差分方法( f d t d ) 的一些研究。全文共分为 四章。 第一章是绪论部分，首先介绍了研究对象和研究背景，特别是时域有限差分 方法的发展历史和发展状况，然后概述了文章的主要内容。 第二章讨论了三维直角坐标系中的f d t d 方法：首先介绍了f d t d 格式以及各分 量节点的取法。然后介绍t f d t d 方法的数值稳定性和数值色散问题：包括时间间 隔的稳定性要求以及c o u r a n t 稳定性条件一一即空间步长和时间间隔之间所需满 足的关系。此外，还介绍了时域有限差分方法的边界条件：先从波动方程角度介 绍了f d t d 的吸收边界条件，由此引出t m u r 的一阶、二阶条件；然后又介绍了一 种吸收效果较好的完全匹配层( p m l ) 吸收边界条件。 第三章从p m l 中采用的离散格式一一指数差分格式出发，在以前的研究基础 上进行了推广，基于t a y l o r 展开定理，使用高阶方法来近似空间偏导数，将指数 差分格式推广到2 n 阶。3 1 节从微分方程基本理论出发，严格推导了时间导数的 指数差分公式。3 2 节根据t a y l o r 展开，将空间一阶偏导数用2 n 个对称点的函数 值的线性组合近似，并根据线性代数理论推导了确定系数的公式。3 3 节采用了 一个有效的变换将指数差分格式转变成一种易于分析的离散格式，并利用 f o u r i e r 方法得到了稳定性条件、阶数与谱域以及阶数与数值色散性的关系。由 于3 2 中推导的2 n 阶格式是一种中心差分格式，对于边界附近的点并不适用，所 以在3 4 节中，讨论了这种情况下的处理方法，并最后给出了数值例子加以验证 算法的有效性。 第四章提出t f d t d 方法的一种新应用：从m a x w e l l 方程和s c r h 6 d i n g e r 方程中 找出关联点一一电流项，然后使用f d t d 格式进行数值迭代耦合求解，从而得到一 种考虑量子效应的电磁场数值计算的新方法。 关键词：时域有限差分( f d t d ) 、完全匹配层( p m l ) 、指数差分、高阶指数差 分、量子效应、耦合数值解。 a b s t r a c t a b s t r a c t t h et h e s i si sm a i n l ya b o u ts o m er e s e a r c ho nf i n i t e d i f f e r e n c et i m e - d o m a i n ( f d t d ) m e t h o df o re l e c t r o m a g n e t i cw a v e t h e r ea r e4c h a p t e r si na 1 1 f i r s t ，w ei n t r o d u c et h eo b j e c ta n db a c k g r o u n do ft h er e s e a r c h ：m a x w e l l e q u a t i o n sa n ds o m en u m e r i c a la l g o r i t h mo nt h e m ，e s p e c i a l l yt h eh i s t o r ya n d d e v e l o p m e n to ff d t o i nc h a p e r2 w ei n t r o d u c et h ef o r m u i a t i o na n dc e l l si nt h e3 - d i m e n s i o nc a r t - e s i a nc o o r d i n a t e s t h e nw ei n v e s t i g a t et h es t a b i l i t ya n dn u m e r i c a ld i s p e r s i o n o ft h em e t h o d 。i n c l u d i n gt h er e q u i r e m e n t sf o rt h et i m ei n t e r v a la n dc o u r a n t c o n d i t i o nf o rs t a b i l i t y , b e s i d e s 。w ef i r s td i s c u s st h ea b o r b i n gb o u n d a r y c o n d i t i o n sf o rf d t df r o mt h ew a v ee q u a t i o n 。t h e nd e r i v et h e1s t - o r d e ra n d 2 r i d - o r d e rm u r sa b s o r b i n gb o u n d a r yc o n d i t i o n ，a tt h ee n do ft h i sc h a p t e rw e p r e s e n ta ne f f i c i e n tb o u n d a r yc o n d i t i o n ：p e r f e c t l ym a t c h e dl a y e r ( p m l ) f o rt h e a b s o r p t i o no fe l e c t r o m a g n e t i cw a v e b a s i n go nt h ed i s c r e t es c h e m eu s e di nt h ep m l v i z e x p o n e n t i a lt i m e d i f f e r e n c e ( e t d ) s c h e m e t h ea u t h o rd o e ss o m ee x t e n s i o na n da d v a n c e m e n t ， e x t e n d i n gt h ec o m m o ne x p o n e n t i a ld i f f e r e n c es c h e m et o2 n o r d e r i nc h a p t e r 3 ，w ep r e c i s e l yd e d u c ee t df r o mt h eb a s i cd i f f e r e n t i a le q u a t i o n st h e o r y t h e n a c c o r d i n gt ot h et a y l o rt h e o r e m ，w eu s et h el i n e a rc o m b i n a t i o n so ft h e f u n c t i o n sv a l u e so n2 ns y m m e t r i c a lp o i n t st oa p p r o x i m a t et h ef i r s to r d e rp a r t i a l d e r i v a t i v e o fs p a t i a lv a r i a b l ea n dg e tt h ed e t e r m i n a t ec o e f f i c i e n t s i n3 3 ， a d o p t i n ga s k i i l f u it r a n s f o r m a t i o nt oc o n v e r tt h ee t oj n 吣as c h e m em a ti se a s y t oa n a l y s ea n dd e r i v et h es t a b i l i t yc o n d i t i o n ，t h er e l a t i o nb e t w e e nt h eo r d e ro f t h es c h e m ea n ds p e c t r u ma sw e l ta st h eo r d e ra n dn u m e r i c a ld i s p e r s i o n a tt h e e n do ft h i sc h a p t e r , w ed i s c u s sh o wt od e a lw i t ht h ep o i n t sw h i c hn e a rt h e b o u n d a r ya n du s ean u m e r i c a le x a m p l et os h o wt h ee f f i c i e n c ya n de f f e c to ft h e h i g ho r d e rm e t h o d ， c h a p t e r4 i sm a i n l ya b o u tac r e a t i v ei d e ao nf d t d w h i c hc o n s i d e r i n gt h e a b s t r a c t q u a n t u me f f e c ti n t ot h ec o m p u t a t i o no ft h ee l e c t r o m a g n e t i cf i e l d t h ep a p e r i n v e s t i g a t et h er e l a t i o n s h i pb e t w e e nm a x w e l le q u a t i o n sa n ds c h r o d i n g e r e q u a t i o n ，a n db r i n gf o r w a r dt h ei d e at oc o u p l e ds o l v et h e mn u m e r i c a l l y k e y w o r d s ：f d t d 、p m l 、e r d 、h i g h - o r d e re t d 、q u a n t u me f f e c t 、c o u p l e d a n dn u m e c i a is o l u t i o n i 1 1 电磁波f d t d 方法的若干研究 第一章绪论 1 8 7 3 年，麦克斯韦( m a x w e l l ) 建立了电磁场的基本方程，微分形式如下 1 ： v x 日：罂+ ， 优 v x ：一娑一，。 o t ( 卜2 ) 其中为电场强度，d 为电通量密度，日为磁场强度，丑为磁通量密度，为电流 密度，厶为磁流密度。在这1 0 0 多年的发展历史中，随着电磁波理论研究与实 际应用的推广和计算机技术的进步，i “a x w e l l 方程的数值解法也得到了大量的关 注和研究。例如矩量法( m o m ) 、有限元法( f e m ) ，边界元法( b e m ) 以及时域有 限差分( f d t d ) 方法 2 ，等等。 1 9 6 6 年，k s y e e 首次提出了一种电磁场数值计算的新方法一一时域有限差 分法( f i n i t ed i f f e r e n c et i m ed o m a i n ，f d t d ) 方法 3 。对电磁场e 、日分量 在空间和时间上采取交替抽样的离散方式，每一个( 或日) 分量周围有四个日 ( 或e ) 分量环绕( 3 维情况下) ，应用这种离散方式，将含时间变量的m a x w e l l 旋度方程转化为一组差分方程并在时间轴上逐步推进地求解空间各点的电磁场 分量。y e e 提出的这种抽样方式就是大家现在熟知的y e e 元胞( y e ec e l l ) ，而 当时他采用的差分格式是一种二阶的显式格式。 在f d t d 的发展历史上，除了y e e 以外，很多数学家和工程学家也作出了大 量杰出贡献。t a y l o r 等用f d t d 分析非均匀介质体的电磁散射问题，提出用吸收 边界来吸收外向行波，吸收边界采用简单差值方法 4 。t a f l o v e 用f d t d 讨论了 时谐场情况的近一远场外推，以及数值稳定性条件 5 。s u i 等提出用f d t d 计算 含有集r f i 参数元件的数字电路模型 6 。b e r e n g e r 提出将m a x w e l l 方程扩展成为 场分量形式，并构成一种全新的吸收边界一一完全匹配层 7 9 。 本义首先介绍了三维直角坐标系中的f d t d 方法的基本形式和取样节点；接 着介绍如j 何从平面波解分析方法的数值稳定性；然后介绍了三种常用的吸收边界 条件：m u r 一阶、二阶吸收边界以及完全匹配层吸收边界；在此基础上引出了指 电磁波f d t d 方法的若干研究 数差分格式，并在第三章中对其进行了推广和改进，得到一种精度更高的高阶指 数差分方法，然后给出数值例子加以验证算法的有效性；此外，这一章还讨论了 高阶方法的谱域、数值色散性与阶数的羌系，稳定性条件，以及边界附近点上的 处理方法等：第四章提出了一种新思想一一如何将f d t d 方法应t f | j 于考虑量子效 应的电磁场数值计算。最后探讨了该方法的研究以及应用的前景和展望。 电磁波f d t d 方法的若干研究 第二章三维直角坐标系中的f d t d 方法 2 1 三维直角坐标系中的f d t d 格式 各向同性介质中的本构关系为二 d = 锺 8 = “h j = o e j 。= 口。h ( 2 一1 ) 这里s 表示介质介电系数，卢表示磁导系数，硪示电导率，而表示磁导率。f 和分别为介质的电损耗和碰损耗。真空中萨o ，= o 以及e = 。：8 8 5 x 1 0 1 2 f m , p 斗。= 4 x 1 0 y m 。 在直角坐标系中，利用上述本构关系将( 卜1 ) 和( 卜2 ) 式写为 ( 2 - 2 ) 警鲁叫警鸣以 誓一等：一卢导吒， ( 2 - 3 ) i 言一卢i 吓乒， 坦。， 誓一等= 一鲁呱丘 对( 2 2 ) 和( 2 3 ) 进行时间和空间的差分离散，令( j ，y ，z ，f ) 代表( e ，e ，e ：) 或 h ( h ，口，：) 的某一分量，记 哦 呜 匾 + + + 峨i 堡西强百 占 占 = = | j 哆i弧i弧砂 幔一钞弧i 鸭i 电磁波f d t d 方法的若干研究 则对时间的差分近似为 f ”( f ，) = f ( i a x ，缈，t 缸，n a t ) ( 2 4 ) 丝学业l 。笪丛学塑( 2 - 5 ) a i ，：m 出 而对空间的差分近似( 以善为例，其余两个类同) 为： 一 。：! ! ! ! ! ：! ：生二：! ! = ! ! ! ! ! ：1 2 血 ( 2 6 ) 图( 2 - 1 ) 就是著名的y e e 元胞( y e ec e l l ) 示意图 3 ，电场和磁场各节点 的空间位置如图所示。从图中可以看到，每个电场分量被四个磁场分量围绕； 同样，每个磁场分量乜被四个电场分量围绕。这种取样方式非常适合f d t d 的 计算，同时也符合法拉第电磁感应定律和安培环路定律的物理结构，并恰当地 描述了电磁波在空间的传播过程。此外，电场和磁场在时间轴上也是交替抽样， 时间间隔为，2 ，这就将i u a x w e l l 方程离散成一组易于迭代求解的显式差分 格式，所以，只要给出初值和边界条件，就可以求出各时刻、各空间位置的电 磁场分布了。 图2 1 f e e 元胞和电磁场各分量节点 一 纠等 电磁波f d t d 方法的若干研究 根据上面的差分格式，我们以尉u 髓为代表写出求解公式 1 0 。( 其余类似) p ”( m ) = c a ( m ) e ；( m ) i 堕型掣蔓尘浯。， 一兰! 兰! ：圭：! ：! ：i ! ：兰! 兰! ! ：i ：! ：! ：墨! 其中肺= ( f + ；，女) ，c a ( ) ： 1 + 口( m ) t 2 f f m l j ( m ) a t 2e ( m ) ：。2 ( m ) = c p ( m ) 日：2 ( m ) 一c q f 竺：! 坐生二竺：尘生。2 - 8 )一( m ) 【 ( ) e 抛，+ 寻，t + 1 ) 一；( f ，+ 导， 其中m ：( f ，j + ；，t + ；) ，c p ( 。) ： 二二 l + 口。( m ) a 2 ( mj 口。( m ) 2 【m ) 若采用上述公式，则可得到时域推进计算流程和计算节点分布 图2 2f d t d 计算流程图 电磁波f d t d 方_ ；去的若干研究 电磁场分量 x y zt 1 + i 2， 膏 f j j + 1 2 膏疗 丘 j k + 1 2 尼 j + 1 2 k + l 2 鞋 足i + 1 2 k 2n + l 2 厦 i 挖 j n 挖 膏 表2 1 计算分量空间、时间分布 值得注意的是，y e e 提出的差分格式，在空间和时问上都采用了中心差分， 根据t a l y o r 展开定理可以得知，这种离散格式的求解精度是二阶的。在后面的 章节中，作者提出了一种高阶的差分格式，在增加计算量的情况下，提高求解精 度到任意2 n 阶。 2 2f d t d 数值稳定性和数值色散 通过2 i 的分析，可以发现f d t d 方法是以一组二阶中心差分方程来代替 m a x w e ll 旋度方程。只有离散后的差分方程组的解是收敛且稳定的，这种数值方 法才是有意义的。本节主要讨论f d t d 方法的稳定性和收敛性对时间和空间离散间 隔的要求。 考虑时谐场情况 1 1 ，即 f ( x ，y ，z ，) = 六e ” 显然这一稳态解是一阶微分方程 蔓：i 研 a 。 ( 2 9 ) ( 2 - 1 0 ) 电磁波f d t d 方法的若干研究 的通解。用中心差分代替( 2 一l o ) 中等式左边的时间微分，得 = c 矿”= j c o f ( x ，y ，z ，n a t )( 2 - 11 ) 当，足够小得时候，定义数值增长因子p 为 则由( 2 - 1 1 ) 和( 2 1 2 ) 可得 ( 2 1 3 ) 的解为 ，”+ j f 一 胪7 5 p 2 一j t o a t p 一1 = 0 p = 竽j 丐 根据数值稳定性要求，当月一佃，出呻0 时。例s l 所以由( 2 1 4 ) 得充分条件 ( 2 1 2 ) ( 2 - 1 3 ) ( 2 1 4 ) ( 2 - 1 5 ) i 由m a x w e l l 方程组可导出，电磁场的任意直角分量均满足齐次波动方程 警+ 等+ 害一专等= 。 缸21 却2 。出2 c 2 西2 一。 根据物理中波矢朋自定义并记t 2 t ，女 ，则有 如e + t ；+ 女；= ( 2 , r t 2 = 等 ( 2 - 1 6 ) ( 2 - 1 7 ) 电磁波f d t d 方法的若干研究 其中拜口，毛分别表示电磁波在介质中的传播速度和波长。所以平面波表达式 ，( x ，y ，：，) = 厂de x p 一j ( k ，x + ，y + i ：z 一删) ( 2 1 8 ) 是( 2 1 6 ) 的解。 用二阶中心差分离散( 2 1 6 ) ，并将( 2 一1 8 ) 代入后得 塑。一s i nz ( - 竽) 影 ( 譬) 2 ( 2 - 1 9 ) 用相似的方法得到( 2 一1 6 ) 中其余两个空间导数的差分格式，则( 2 - 1 6 ) 的离散近似 可以写成： 警+ 警十等等。 。， ( 等) 2( 等) 2( 等) 2 ， 移项，然后两边同乘以( 孚) ：得到 竺兰+ 竺兰+ c 争2c 爷 ( 警 由条件( 2 1 5 ) 和不等式放缩得到c o u r a n t 稳定性条件( 这是一个充分条件) ： c a t 2 它给出了空间网格和时间步长应该满足的关系。 特别的，当血= 缈= 止= j 时，有 ( 2 - 2 2 ) = 1【- 电磁波f d t d 方法的若干研究 a 在三维情况下，f f 牟，即时间步长应小于等于波以光速通过y e e 元胞1 3 3 对角线的所需时间， b 在二维情况下，c 址s 皂，即时间步长应小于等于波以光速通过y e e 元胞1 2 2 对角线的所需时间。 c 在一维晴况下，c 址s 占即时间步长应小于等于波以光速通过y e e 元胞的所 需时间。 此外，即便介质本身无色散，但对波动方程做差分近似的做法本身会引入波 的色散，这种现象称为数值色散。而且使用f d t d 计算，会出现波的传播速度与 传播方向相关，这也是由数值离散引起的，称之为各向异性。虽然这些数值方法 本身的误差无法彻底逶免，但是可以通过各种分析手段得到尽量小的误差。一般 情况下要求( 详细推导见文献【1 2 ) ： 缸，妙，血詈，出吾 ( 2 2 3 ) 这里旯和7 分别是电磁波在介质中的波长和周期。 通过以上的分析，我们了解到要保证f d t d 方法的稳定性和收敛性，对时间和 空间的取样间隔都有一定要求，总结起来主要有以下三点： 1 ) 时间步长的选择应该满足c o u r a n t 条件，以保证数值稳定性。 2 ) 空间间隔和时间步长的选择应尽量减小数值色散。 3 ) 空间间隔应满足一定条件以减小各向异性。 2 3f d t d 的边界条件 由于计算量和存储量的限制，f d t d 计算只能在有限区域进行。为了能模拟 开域中的电磁传播过程，在计算区域的截断边界处必须给出吸收边界条件。关于 f d t d 吸收边界的工作也是该方向的一个研究热门，从开始简单的插值边界【1 3 】， 到后来广泛采用的m u r 吸收边界 1 4 】，以至9 0 年代中后期的完全匹配层( p m l ) ， 电磁波f d t d 方法的若干研究 其吸收效果越来越好。此外，为了适应不同介质和应用模型各种相应的高性能 吸收边界也层出不穷 1 5 1 8 】。在这章中我们就最常用的边界条件进行了分析 和讨论：包括m u r 的一阶、二阶吸收边界和p 扎吸收边界。特别是最后介绍的 p m l 吸收边界中使用的一种差分格式是第三章内容的背景和基础，而本文的一个 工作就是对这种差分格式进行了改进和推广，具体内容详见第三章。 从m a x w e l l 方程组出发，可以导出电磁波的电场f 和磁场所满足的波动方 程，在二维情况下其表达式为： 在无界均匀媒质条件下，( 2 2 4 ) 的解为平面波解 1 1 ： 其中 ( 2 2 4 ) f ( x ，y ，t ) = a e x p j ( 耐一kx 一女j ，) ( 2 2 5 ) 女：+ 女：= k2 = ( 竺) z 。 c 将( 2 - 2 5 ) 代入( 2 - 2 4 ) 并保留对x 的二阶导数，得： 记 粤+ ( k 2 - k ；) f ：o 咖。 卜昙一，瓜 k 昙+ ，厩 则有：( 工+ 工一) ，= 0 ( 三一为左行波算子，0 右行波算子) 。 ( 2 2 6 ) ( 2 - 2 7 ) ( 2 2 8 ) 划 塑酽一 ， 塑矿 + 盟掰 电磁波f d t d 方法的苦干研究 设z = x 。平面为截断边界，且在计算区域的左侧。因为在计算区域中同时存在入 射波和反射波： 厂一e x p + i 二再+ 删 ( 2 2 9 ) ，。= a + e x p j ( 耐一厣飞x + k 帆( 2 - 3 0 ) 这里正为左行波，即入射波：而工为右行波，即反射波，且满足，= 工+ 正。 若在截断边界处设置条件：一一。= o ，就相当使反射波成分为o ，再根据频域 到时域的算子转换，得到左边界上的吸收边界条件 1 2 】： d l - q - - 一 “ 同理可得右边界上的吸收边界 c 晏+ ( 2 3 1 ) ( 2 3 2 ) 对上述的吸收边界条件相应的频域形式使用t a y l o r 一阶、二阶级数展开，再转 换到时域表达式就得到了m u r 一阶、二阶吸收边界条件( 以左边界为例) 1 2 】： 一阶： 二阶： c 去一与t = 。 s 。， c ；岳专箬+ 三蒡叫。= 。 s a ， ：丽一7 矿+ j 萨，l 划 ( 2 书4 ) 完全匹配层( p e r f e c t l ym a t c h e dl a y e r ，p m l ) 首先由b e r e n g e r 于1 9 9 4 年 提出 1 9 。通过在f d t d 区域截断边界处设置一种特殊的介质层，该层介质的波 阻抗与相邻介质波阻抗完全匹配，因而入射波将无反射地穿过分界面而进入p m l 电磁波f d t d 方法的若干研究 层。实际计算证明，p m l 对入射波有很好的吸收效果 2 0 】。 在二维t e 波的p m l 介质层中，b e r n g e r 假设将磁场分量：分裂为两个子分 量h 。和h 。，且有：= h ，+ h ，进而l a x w e l l 方程转化为： s 鲁坞t = 坚裂 占拿u e ，： d f 。 a ( 日。+ h 。) 孥+ 盯。以： 。 敏 岫。 “等岭 口v 葩。 出 0 e 卯 ( 2 - 3 5 ) 其中盯，盯，盯。，盯，为介质的电导率和导磁率。定义介质的波阻抗为电场和磁场 之比； ，e1 z 2 百2 i 石 ( 2 3 6 ) e ，u 分别表示p m l 介质的介电系数和磁导系数，则为了达到吸收效果，这些参数 必须满足p m l 的重要基本条件一一阻抗匹配条件： f 生：垒 l e 1 生：垒 【 ( 2 - 3 7 ) p m l 吸收层的网格剖分方式与常规的f d t d 网格完全致，唯一的区别是由于吸 收介质中的波方程将h = 分量分裂为日。和h 。，因而需在原来也分量的节点处 同时计算- ，和h 。由于在p m l 介质中电磁波衰减很快，常规的f d t d 中y e e 的 差分格式已不再适用【1 9 以2 】，方程( 2 - 3 5 ) 的离散将采用指数差分 2 0 ： 电磁波f d t d 方法的若干研究 掣“+ ；) = e x p ( - q ( ，+ 圭) 譬) 硪“+ 圭) 1 一e x p ( 一盯，( ，+ 芝1 ) a 占t ) 一( + 吉) n + - l n + - i 型坠掣】 ( j + j 1 ，) = e x p ( 一吒( f + 互1 ) a s t ) 啪+ j 1 ，) 一譬，。绁n,t_等啦11 ， q ( f + = ) 血 1 h ：；( f ，) ：e x p f 一盯。( f ) 竺) 日，( f ，- ，) 一荨a t1 竺趱5 x ， 盯眦( f ) 。 ) _ e x p ( - 叫砖) 九，) 一! ：! ! ：! 竺差- 竺! ! ：! ：皇二竺竺! 二兰一 些r ： 2 ：! 仃。，( f ) 。 缈 特别的，对v 口，= o ，j = m x ，m y ，工，y ，：t - j 血或缈有 1 1 - e x p ( - o a t ) i f zm1 上= 兰：兰 q ”d i o i s 2 1 - e x p ( 嵋丝) “fl a u2 z ( 2 3 8 ) ( 2 - 3 9 ) 这里，我们不加推导地给出了p m l 介质中的指数差分格式，这其实是一种二阶指 数差分。在下一章内容中，严格地推导了指数差分的公式，并将其改进、推广到 了2 n 阶，得到高阶指数差分的一般形式。 电磁波f d t d 方法的若干研究 3 1 引言 第三章高阶指数差分方法 前一章讲到的三种吸收边界条件中，效果最好的是完全匹配层( p 札) 吸收 边界。b e r n g e r 在他提出p m l 方法的文章 1 9 中指出，对于损耗介质，应采用不同 的计算方法以达到更好的近似效果，因而他提出在p m l 吸收层中使用对时间导数 的指数差分格式，在文献 2 0 2 2 中也对这种离散格式进行了分析和讨论。然而 大部分支献中一般只是在p m l 中采用这种格式，而且使用的都是二阶指数差分， 并采用场量分裂方式( 根据p m l 的特征) ，若这样直接应用到整个计算区域，则不 仅对存储量要求高，而且精度低。文献 2 3 推导了四阶指数差分公式，并将其使 用在损耗介质的整个计算区域中，这种差分格式在不增加存储需求的前提下，提 高了计算精度，达到了更好的数值效果。 本文在这些研究基础上进行了一些改进和推广，基于t a y o r 展开定理，使用 高阶方法来近似空间偏导数，将指数差分格式推广到任意2 n 阶。在这一章中，首 先从微分方程基本理论出发，严格推导了时间导数的指数差分公式。在第二节中， 根据t a y l o r 展开，将空间一阶偏导数用2 n 个对称点的函数值的线性组合近似，并 从线性代数理论中推导了确定系数的公式。第三节，作者使用了一个巧妙的变换 将指数差分格式转变成一种易于分析的离散格式，并利用f o u r i e r 方法得到了稳 定性条件、阶数与谱域以及阶数与数值色散性的关系。由于第二节中推导的2 n 阶格式是一种中心差分格式，对于边界附近的点并不适用，所以在第四节中，作 者讨论了这种情况下的处理方法。在本章的最后给出了数值例子加以验证算法的 效果。 电磁波f d t d 方法的若干研究 3 2 指数差分 在大部分文献 1 9 2 2 中，指数差分仅作为一种离散格式直接给出，为了更 好地理解这种方法，在这一节中作者给出了严格的数学推导。以( 2 2 ) 中第一 个式子为例，我们将它重新整理为： 堡+ 吐o - ；上掣一三掣( 3 _ 1 ) o t3 却e 跣 记g ( ) ：望熹一鬈则根据微分方程理论，( 3 - 1 ) 的解t ( ，) ：e x 以) + 瓦( ，) 洲口 o - do - 其中基本解量。( f ) = c e ，特解三。( f ) = 0k ( f 弦j f 出弦j ，c 为任意常数。所 以得到： 记e ：= t ( n a t ) ，贝0 删= ( 1 l g ( ) e ：7 出+ c ) e ( 3 2 ) ( 3 3 ) e = ( i 1j g c r ，e e r t 衍+ c 4 ，：。e “ c 。一a ， 一! ， 由( 3 - 3 ) 一( 3 4 ) e5 可得： f “母e 一拶4 舭却“ s ) 出 肿 f f 卅 = o + 出b i 传 0 f = + e 电磁波f d t d 方法的若干研究 所以只要对定积分 ：“”g ( r 弦_ ，d t 使用不同的数值积分方法就可以得到不同精 叫 度的解，例如： a 直接取中间值 b 梯形公式 r 徘g ”+ 蛩。 ( 。_ 6 ) 。n + 1 ) a t 鼎弋半。憎e 泞7 ) 鹋咖湖公式f 净腓。钞。竿+ 。掣“吲。删 其中( 3 - 6 ) 就是常用的指数差分格式的一般形式。 3 3 高阶空间差分离散 在3 2 中我们已经推导了指数差分方法的一般形式，如引言中所述，在微 分方程的另一边，我们可以基于t a y l o r 展开，将电磁场分量对空间的导数用高阶 t a y l o r 差分格式去逼近以达到更好的近似效果。 定理3 1 ：如果函数，( x ) 满足条件：在点j = 的某邻域u ( x 。，d 内有定义； 在此邻域内存在n + l 阶导数；在z = 处有n 阶导数，那么，( z ) 在u ( x 。，国内 可以写成下面的形式： 厂( x ) = ；| ；。：! ：i 屿( i ( 工。) + 。( ( x x 。) ”1 ) ( 3 9 ) 这就是著名的t a y l o r 级数展开定理，利用这个定理我们来推导具有2 n 阶精度的 差分格式。 以善方向的偏导数为例，f ( x ) 在x = ( ”+ f ) 缸= 至字，册= l ，2 ，) 各 点上取值，共2 n 个点。将这些点的函数值在z = i a x 处t a y l o r 展开，则可得到 电磁波f d t d 方法的若干研究 矩阵表达式 ( 3 - 1 0 ) 其中： f = 彳口+ d ( ( x ) 2 ”“) f = 听。正。，；：。) 。厶：。) 2 】0 。 口= 【兀，矗”，嚣“，矗”，砖2 川如州。 4 川2 + 1 1 = 1 1 a t 2 一x 2 3 a x 2 3 z 2 ( 血2 ) 2 2 1 ( - 缸2 ) 2 2 1 ( 3 z k 2 ) 2 2 1 ( - 3 缸2 ) 21 2 1 1 ( 2 n 1 ) z k 2 ( 2 n 一1 ) 5 x 2 2 1 1 - ( 2 n - 1 ) a x 2 - ( 2 一1 ) 四2 2 1 ( 缸2 ) 2 ”( 2 | v ) ! ( - 缸2 ) 2 ”( 2 j ) v ) ! ( 3 6 x 1 2 ) 2 n ( 2 叼! ( - 5 x 2 ) 2 ”f 2 忉! ( 2 n 1 ) a x 2 “( 2 忉! 【韶n l 刎2 ”( 2 n ) ( 3 1 1 ) ( 3 一1 2 ) ( 3 1 3 ) 口f ( 缸) 2 肌1 ) = 【d ( ( 缸) 2 肌1 ) ，o ( ( 血) 2 州) ，o ( ( 血) 2 川) 。 ( 3 1 4 ) 因为在m a x w e l l 方程中仅涉及电磁场分量的一阶偏导数，所以为了避免计算高阶 导数，就可利用多个节点函数值的线性组合来表示要l ，以达到任意阶精度。 c i 卢，缸 而各个系数可以用求解线性方程组的办法来得到。 要l = 去。赘伽叫以叫一。c c 则由( 3 - 1 0 ) ( 3 - 1 5 ) 可得到系数确定的方法，由于采用了中心差商，所以令 d ( z 。叫，：= 一以枷，：，则所有偶次项导数都可消去，为求得系数我们可以求解以 下齐次线性方程组： 电磁波f d t d 方法的若干研究 ( 1 2 ) 3 ( 3 2 ) 3 【( z 一1 ) 2 3 0 2 ) 5 ( 3 2 ) 5 掣一0 2 5 j ! l 2 ) 2 ”4 ( 3 2 ) 洲( c 2 一1 ) 2 】2 ( 3 一1 6 ) 得到0 1c 2 o 一。) 后( 显然系数矩阵的行列式不为0 ，故有唯一解) ，经整 理得： d t z ，一- ，z = 一d t z 。一- ，z = 三j i = j 彳x ! ；：三；：；：；= 三： ；i ；j i i 二i f c s 一， 这一结果与文献 2 4 2 5 是一致的。 3 4 高阶指数差分方法 根据3 2 和3 3 的分析，我们将一般的指数差分格式进行推广一一得到高 阶指数差分方法。具体地说就是对时间偏导数使用指数差分格式，而对空间偏导 数使用高阶中心著分展开，从而得到相对较好的数值效果。具体公式如下( 以( 3 1 ) 为例) ： 掣l ( f 十j 1 川= e 珊+ 如卅 + 季专羔，吼。州：确+ 抄莩问， 一古。姜+ m ：删，：日，c ，+ j 1j 七+ 莩) 】) ( 3 - 1 8 ) 其余几个分量的公式类似。 rl，k 电磁波f d t d 方法的若干研究 首先，我们以一维的高阶指数差分为例，进行稳定性条件讨论( 因为所做的 处理都是针对时间变量的，所以可以容易地推广到二维或三维) ： 蹦一吲力一譬鸯m ”掣仔 硝岭i 硝呤芸。赫珊圳伊z 这里表示空间网格大小，= ，- 。= 孑瓦币岽三蔷蕃 为了易于对这种高阶指数方法进行分析，我们首先引入一个变换： 变换( 3 - 2 1 ) 将z 【o ，m ) 寸z 【0 , 2 】，记 ( 3 2 1 ) f ：三_ f - ： 盯盯 f 。：旦寸f 。：卫 d m盯m ( 3 2 2 ) ：坐_ r ：些 ff k ：坐一，k - ：g 则利用( 3 - 2 1 ) 、( 3 - 2 2 ) 将( 3 - 1 9 ) 和( 3 2 0 ) 变为： h a t 础轳矗暇素珐莹叱一m 蛳n + ；2 m ：- 1 净z 。， 蔓 = 电磁波f d t d 方法的若干研究 ，h 。竺 确”扣差。一哆惫，t 去。弘刚m z a ， 至此高阶指数差分经变换从形式上与中心差分格式一致。在2 章中，我 们介绍了应用平面波解得到的f d t d 方法的稳定性要求，这里我们将使用微分方程 数值解中常用的f o u r i e r 方法 2 6 来进行高阶指数差分算法的稳定性分析： ff ( f ) 1 令巩2 l h ( f + 争j 则巩“2 g 面月，其中g 为一个2 2 矩阵，各个元素为 瓯= 等一筹萎n 驰n 。k 。s i n ( 丁2 m - 1 抛胁( 莩眦) 】 g l z 一焉酶n 。ms i n ( 丁2 m - i 抛， g z ，一，2 ，i a t 争篙d l z i ) 2s i n t 2 m - 1 g 2 22 笋 这影班班钞= 等= 畚_ 表示虢 经分析得稳定性条件： 再利用变换( 3 - 2 1 ) ，得到2 n 阶指数差分方法的稳定性条件为： ( 3 - 2 5 ) ( 3 2 6 ) 兰c 兰多 c s z ， 这里霄表示电磁波所在空间维数，所以胆1 ，2 ，3 。特别的，设= 1 且h 专0 ， 则因为嬲；( 兰苦) = 1 ，萎nk ：。m f = l ，所以此时( 3 2 7 ) 退化成c o u r a n t 稳定性 亡。 上压 一 p 南 上压 一 电磁波f d t d 方法的若干研究 条件。这说明( 3 - 2 7 ) 是差分方法稳定性条件的一般表达式，它包含c o u r a n t 稳定 性条件这种特殊情况。 同样利用f o u r i e r 方法( 又称v o nn e u m a n n 方法) ： 厂( 0 + 提) ) = f ( s a ) e x p “6 得到高阶指数差分方法的谱域表达式 3 2 8 ) 拈丢乳。珧i n ( 2 m 厂- 1 ) 】 ( 3 _ z 9 ) 即用数值频率云来近似理想频率k 。图3 1 给出了2 、4 、6 、3 2 阶等几种情况下的 谱域特性比较，可以看出，指数差分阶数越高，在更宽的谱域内对理想频率的逼 近程度就越好，因而得到的计算结果误差就越小。 图3 1 不同阶数指数差分的谱域特性 利用损耗介质参数的复介电常数处理，可以得到高阶指数差分色散关系的表 达式。事实上，根据文献1 2 7 ，对损耗介质定义：尹= 占l ，昙，则这种离散格式 在一维情况下的数值色散关系为： 电磁波f d t d 方法的若干研究 兰【吐in【2m，-i云2m-12s l m - - l 云- ) 】【吐i n 【， ) 】 2 对于二维有 t = f c o s 8 ，t = k 一 s i n o ( 3 - 3 0 ) ( 3 - 3 1 ) 对于三维有 丘t s i n g c o s # ，t = k 一 s i n 8 s i n ，t = 云c o s 8 ( 3 3 2 ) 图3 2 给出了二维情况下不同阶数指数差分方法的数值色散特性一各向异性， 从图中可以看出阶数的提高可以有效减小数值色散的各向异性，但超过6 阶时变 化已经不大。 o 们 i0 0 】 1 伽o 氆9 9 3 od 9 争d 9 9 辞 0 9 9 2 0 9 9 0 09 聃 o 朔6 0 9 01 8 0 2 7 0 3 6 0 嚣，( ) 图3 2 不同阶数指数差分的数值色散特性 ( 血= 缈= ，丑= 1 0 ，一= 1 歹万，矿表示波的相速度，目为传播方向角) 从公式( 3 - 1 8 ) 可以发现，要使用2 n 阶方法，相关分量的空间下标必须满足 l i 2 n 2 - 1 ，f + 型笋】d ，这里。表示整个2 n 阶指数差分方法的计算区域，所 ，p譬： 电磁波f d t d 方法的若干研究 以对v m 【一+ 1 ， & ，一竺冬兰d 即边界附近的点，原方法都会失效。为了 二 在这些点上保持同样的计算精度，我们使用向前、向后差分方法，也就是用 ( 2 n + 1 ) 个计算区域内的点的函数值的线性组合来近似一阶空间导数，使得近似 阶数依然是2 n 阶。由于推导方法类似，在此就仅给出计算公式。详细摧导过程见 文献 2 4 2 5 。 对于向前差分，我们有： 其中： ：_ 1 2 ng 。a + ，+ o ( a 2 ”) 凸z 0 ( 3 - 3 3 ) f 一i t = o 舻rc s 嘎， l ! 二! 堡丝： 1 2 阶指数差分。 ( 4 ) 某些电磁波的应用、计算和模拟中，对运算精度的要求非常高，因而需采用 高阶方法 2 8 。 综上所述，对空间导数采用高阶的差分格式还是有一定意义的。 电磁波f d t d 方法的若干研究 4 1 引言 第四章考虑量子效应的f d t d 方法 近年来，材料科学、半导体以及电子产业技术迅猛发展，产品尺度已到纳 米级水平。这就需要使用新的方法来分析电子器件及集成电路的各种新特性一一 即将量子效应考虑于其中。例如，在量子点( q d ：q u a