（通信与信息系统专业论文）高阶精度时域有限差分方法的研究及其应用.pdf

摘要 摘要 时域有限差分方法有着良好的瞬态性和广泛的适用性，已被广泛应用于各类 电磁场问题的求解中。传统时域有限差分方法在时域和空域上都采用2 阶精度的 中心差分方法，其计算精度较低，数值色散误差较大。对电大问题作电磁波传播 长期响应分析时，由于误差的积累，往往造成波形的严重失真。这是传统时域有限 差分方法的固有缺陷。 为改进时域有限差分方法的精度特性，近年来学者们提出了许多高阶精度的 时域有限差分方法。其中，基于辛策略的高阶时域方法有着优良的精度特性。本 文采用空域上的离散奇异卷积方法，通过改变单边计算带宽m 的值，使得空域精 度达到2 m 阶。时域上采用辛积分传播子方法，其精度可达4 阶。针对时域有限差 分方法的计算精度和计算效率，做了以下创新性和探索性研究。 ( 1 ) 推导了基于波动方程的高阶算法。该方法利用离散奇异卷积方法对波动 方程的2 阶偏微分进行差分近似，空域精度可达2 m 阶。同时利用辛积分传播子方 法构造4 阶时域精度。利用数值算例分析了波动方程形式的高阶算法的精度特性， 同时分析了该算法在计算复杂边界问题时的灵活性。以规则波导的解析解作为参 照，本论文的高阶算法比传统的时域有限差分方法的误差低近2 0 0 d b 。同时计算了 脊波导的t m 模式下的前几个截止波数，该结果与已有文献的结果吻合。 ( 2 ) 推导了基于辛策略高阶时域有限差分方法的完全匹配层的表示形式，同 时针对该算法的特点修改了总场散射场体系的连接边界条件。基于以上结论，利 用基于辛策略的高阶时域有限差分方法计算目标的电磁散射问题。计算了单个导 体方柱、多个导体柱的双站雷达散射截面。数值结果表明，在计算目标的散射特 性时，高阶算法与传统时域有限差分方法的计算结果吻合，但高阶算法使用的网 格数低于普通时域有限差分方法，从而有着更高的计算效率。同时，利用周期边 界条件，计算了一维金属带周期阵列结构的单站和双站雷达散射截面，验证了该 结构的吸收电磁波的特性。同时，进一步计算了覆盖有一维金属带周期阵列结构 的金属导体柱的散射特性。结果表明，覆盖一维金属带周期阵列结构后前向散射 有所增强，而后向散射得到衰减。 摘要 关键词：时域有限差分，离散奇异卷积，辛积分传播子，完全匹配层，雷达散射 截面 i i a b s t r a c t 一一一_ a b s t r a c t t h ef i n i t e d i f f e r e n c et i m e - d o m a i nh a sb e e nw i d e l yu s e di ns o l v i n g e l e c t r o m a g n e t i c p r o b l e m sb e c a u s eo fi t sg o o dp e r f o r m a n c eo fi n s t a n t a n e i t ya n da d a p t a b i l i t vw h e ni t c o m e st ot h et r a d i t i o n a lf i n i t e d i f f e r e n c et i m e d o m a i nm e t h o d ，s e c o n d o r d e rc e n t r a l d i f f e r e n c es c h e m e ，w h i c ha r i s e sl o wa c c u r a c ya n dg r e a td i s p e r s i o ne r r o r , i si n v o l v e di n b o t ht h es p a c ed o m a i na n dt h et i m ed o m a i n f o rt h el a t e t i m ea n a l y s i so f e l e c t r i c a l l v l a r g ed o m a i n s ，t h ew a v e f o r mi so f t e nd i s t o r t e ds e r i o u s l yb e c a u s eo ft h ea c c u m u l a t i o no f n u m e r i c a le r r o r s t h i si st h ei n h e r e n td r a w b a c k so ff i n i t e - d i f f e r e n c et i m e d o n l a j n m e t h o d t oi m p r o v et h ea c c u r a c yp e r f o r m a n c eo ff d t dm e t h o d ，m a n yh i g h o r d e rf i n i t e d i f f e r e n c et i m ed o m a i nm e t h o d sh a sb e e np r e s e n t e dd u r i n gt h e s ed e c a d e s a m o n gt h e s e h i g ho r d e rs c h e m e s ，t h et i m ed o m a i ns c h e m eb a s e do ns y m p l e c t i cm e t h o dh a se x c e l i e n t p e r f o r m a n c eo fa c c u r a c y i nt h i sw o r k ，t h ed i s c r e t es i n g u l a rc o n v o l u t i o nm e t h o di s i n v o l v e di na p p r o x i m a t i n gt h es p a t i a ld e r i v a t i v e sw h e r e 删ho r d e ra c c u r a c yc a nb e a c q u i r e d m e a n w h i l e ，t h es y m p l e c t i ci n t e g r a t o rp r o p a g a t o rm e t h o di sa p p l i e di nt i m e d o m a i nw h e r et h ea c c u r a c yi nt i m ed o m a i ni si m p r o v e dt o4 t ho r d e nt ot h ei s s u eo ft h e a c c u r a c ya n de f f e c t i v e n e s so ff d t d ，i n n o v a t i v ea n de x p l o r a t o r yw o r ka r el i s t e da s f o l l o w 1 h i g ho r d e rf d t dm e t h o db a s e do nw a v ee q u a t i o n si sd e r i v e d 。w h e r ed i s c r e t e s i n g u l a rc o n v o l u t i o nm e t h o dw i t h2 m - t ho r d e ra c c u r a c yi si n v o l v e di na p p r o x i m a t i n g t h e2 n 。o r d e rp a r t i a ld e r i v a t i v e s a n dt h es y m p l e c t i c i n t e g r a t o rp r o p a g a t o rm e t h o di s a p p l i e di nt i m ed o m a i nw h e r e4 t ho r d e ra c c u r a c yi sa c h i e v e d n u m e r i c a le x a m p l e ss h o w t h ea c c u r a c yp e r f o r m a n c eo ft h ep r e s e n t e d h i g ho r d e rs c h e m e ，a n dt h ef l e x i b i l i t yi sa l s o p r o v e di na n a l y s i so fc o m p l e xb o u n d a r yc o n d i t i o n s r e f e r r i n gt ot h ee x a c ts o l u t i o n so f r e g u l a rw a v e g u i d e s ，t h ee r r o r so fh ig ho r d e rs c h e m e si sa b o u t2 0 0 d bl e s st h a i lt h e t r a d i t i o n a lf d t dm e t h o d m e a n w h i l e ，t h ef i r s t4c u t o f fw a v en u m b e r so fr i d g e d w a v e g u i d e sa r ea c q u i r e db yo u rh i e , ho r d e rs c h e m e s ，w h i c hc o i n c i d ew e ! w i t ht h e r e s u l t sa l r e a d yb e e np u b l i s h e d i i i a b s t r a c t 2 t h ea p p l i c a t i o no fp e r f e c t l ym a t c h e dl a y e rt o ( 2 t4 ) s c h e m ef d t di sd e r i v e d ， w h e r e2 m - t ho r d e rd i s c r e t e s i n g u l a rc o n v o l u t i o ni ns p a c ed o m a i na n d4 h o r d e r s y m p l e c t i ci n t e g r a t o rp r o p a g a t o rm e t h o di nt i m ed o m a i ni sa p p l i e d m e a n w h i l e ，t h e c o n n e c t i n gb o u n d a r yc o n d i t i o no ft o t a l f i e l da n ds c a t t e r e df i e l db a s e do nh i g ho r d e r s c h e m ei sa l s op r e s e n t e d b a s e do nt h ec o n c l u s i o na b o v e ，t h eb i s t a t i cr c so fb o t h s i n g l ec o n d u c t i n gc y l i n d e ra n dd o u b l ec o n d u c t i n gc y l i n d e r sa r ec o m p u t e d t h er e s u l t s s h o wt h eg o o dp e r f o r m a n c eo f ( 2 m , 4 ) s c h e m ef d t d a p p l y i n gt h ep e r i o d i cb o u n d a r y c o n d i t i o n s ，b o t ht h em o n o s t a t i ca n db i s t a t i cr c so fo n e d i m e n s i o n a lm e t a l l i cw i r ea r r a y a r ec o m p u t e d ，a n di t sp e r f o r m a n c eo f a b s o r b i n ge mw a v e si sp r o v e d m e a n w h i l e ，t h e s c a t t e r i n gp r o p e r t i e so fm e t a l l i cc y l i n d e rc o a t e db ym e t a l l i cw i r ea r r a ya r ea n a l y z e d r e s u l t ss h o wt h a tt h ef o r w a r ds c a t t e r i n gi ss t r e n g t h e n e dw h i l et h eb a c k w a r ds c a t t e r i n g i sd e c r e a s e dt os o m el e v e l k e y w o r d s ：f i n i t e - d i f f e r e n c et i m e - d o m a i n ，d i s c r e t es i n g u l a rc o n v o l u t i o n ，s y m p l e c t i c i n t e g r a t o rp r o p a g a t o r , p e r f e c t l ym a t c h e dl a y e r , r a d a rc r o s ss e c t i o n i v 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作 及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为 获得电子科技大学或其它教育机构的学位或证书而使用过的材料。与 我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的 说明并表示谢意。 签名： 多、丹， 关于论文使用授权的说明 本学位论文作者完全了解电子科技大学有关保留、使用学位论文 的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁盘， 允许论文被查阅和借阅。本人授权电子科技大学可以将学位论文的全 部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描 等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后应遵守此规定) 签名：壬农9 ，导师签名：盔丝乏主 日期：狮呓年6 月6e t 第一章引言 第一章引言 1 1时域有限差分方法的研究进展 2 0 世纪6 0 年代以前，电磁场理论和工程中的许多问题大多采用解析或渐进的 方法进行处理，这种方法能够得到问题的准确解，而且计算效率比较高，但适用 范围较窄，只能求解具有规则边界的简单问题，对任意形状的边界则无能为力或 需要很高的数学技巧。2 0 世纪6 0 年代以后由于电子计算机的迅速发展，使大型数 值计算成为可能。自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 e d o m a i n ，f d t d ) 方法以来，这一方法很快就引起了人们的关注并得到了很快的发 展。f d t d 的优越性能，结构简单，使它广泛的应用于电磁散射、电磁兼容、微波 电路以及各类导波系统计算问题等电磁场工程的各个方面中。时域有限差分方法 的发展如下所述。 1 9 6 6 年，y e e 提出了直接在时域中求解m a x w e l l 方程的f d t d 方法【l 】。1 9 7 5 年，t a f l o v e 等提出y e e 算法的数值稳定条俐2 1 。1 9 8 1 年m u r 提出了基于y e e 网格 的二阶吸收边界条件【3 1 。1 9 8 5 年t a f l o v e 等提出了f d t d 计算电磁散射的模型【4 1 。 1 9 9 0 年建立了微带线的f d t d 模型【5 】。1 9 9 2 年建立集总参数元件的f d t d 模型【6 】。 19 9 4 年b e r e n g e r 提出了二维f d t d 网格的完全匹配层( p e r f e c t l ym a t c h e dl a y e r s ， p m l ) 呖j t 收边界条件【_ 7 1 ，之后又将二维p m l 吸收边界条件推广至三维计算中【8 j 。1 9 9 6 年柳清伙等人提出时域伪谱方法( p s t d ) 【9 】，同年k r u m p h o l z 提出基于小波方法 的时域多分辨率( m r t d ) 方法【1 0 1 。1 9 9 9 年提出了无条件稳定的交变隐式差分方 向的f d t d 方法【l l 】，之后并将该方法推广至三维问题【1 2 】。2 0 0 0 年以后，f d t d 方 法与其他方法的结合日益增多，并行算法也日益被重视。 另一方面，尽管时域有限差分方法应用日益广泛，但传统的f d t d 法由于受 到数值色散条件对空间网格尺寸的限制，其精度特性无法与基于格林函数的积分 方法相比。高阶f d t d 方法在解决这一问题方面具有较广阔的前景。高阶f d t d 包括空间的高阶差分和时域上的高阶差分。空间高阶差分大多是采用采样函数构 造高阶插值方程得到的。其中比较常用的有如( 4 ，2 ) 精度的f d t d 1 3 】，其基本 原理是通过拉格朗日多项式展开后获得高阶差分方程，在空间上可以得到4 阶精 度。邵 1 4 】【1 5 】等提出了一种通用的高阶空间插值方式，即离散奇异卷积方法。该方 电子科技人学硕士学位论文 法基于离散奇异卷积( d i s c r e t es i n g u l a rc o n v o l u t i o n d s c ) 的d e l t a 核对空间的场量 进行采样并构造高阶差分格式，可以达到很高的精度。d e l t a 核的种类很多，可以 是s h a n n o n 核，p o i s s o n 核以及拉格朗日核。实际上使用拉格朗日核的离散奇异卷 积方法就是使用拉格朗日多项式插值。另一种常用的高阶算法是多尺度分辨率的 时域方法( m r t d ) 如前所述。该方法基于b a t t l e l e m a r i e 紧尺度函数将电场磁场 分量展开，用h a r r 小波尺度函数对时间展开。这些高阶算法对数值色散都有很好 的抑制作用。 同时，在时域上，t a k u oh i r o n o 【1 6 】【1 7 】等提出了一种使用辛积分传播子方法构造 的4 阶时域精度的方法。该方法用辛积分构造显示辛格式，从而构造不同阶数的 辛算法。最早提出辛算法的是冯康先生【l8 1 ，他于1 9 8 4 年国际几何和微分方程会议 上系统地提出了一种能够保持h a m i l t o n 系统的基本特征的算法辛算法。近年 来关于辛算法的文献详见 1 7 j9 1 。另一种时域高阶算法是基于龙格库塔方法的显示 时域方法 2 0 】。该方法使用龙格一库塔方法对第n 步至第n + l 步进行差分离散，构造 高阶差分公式。文献【2 1 】以矩形谐振腔为算例，比较了不同阶数的辛算法和龙格库 塔方法的精度特性。文献【2 2 】讨论了不同阶数的辛算法的数值色散特性。 变步长网格的研究具有重要的意义。对于电大、电小尺寸问题，均匀网格己 不再适用，因为为了保证局部细微结构的精度，网格的大小受到限制，过多的网 格导致内存开销过大，网格尺寸过大则导致精度的降低。因此，变步长网格已成 为迫切的需要。对于变步长网格，较为常见的有两大类，一类是亚网格技术，即 计算区域分为若干块，细微结构处使用细网格，场量变化缓慢处使用粗网格，不 同块之间进行信息交换【2 3 】【2 4 】。另一类是均匀渐变网格，网格之间的尺寸按一定的 规律分布 2 5 。2 7 1 。对于非均匀网格技术，其应用前景较为广泛。尤其是对于电大尺 寸目标的散射问题、超细微结构的微波电路的计算问题都有着很好的计算效率。 由于传统的f d t d 方法受c o u r a n t 稳定条件限制，使得f d t d 的应用范围受 到限制。为了保证细微结构的精度，空间步长必须很小，使得时间步长也相应变 小，计算时间猛增。为此，t n a m i k i 等将交变隐式差分方向方法应用于f d t d 方 法，提出( a d x f d t d ) 方法，将其应用于二维t e 波问题的模拟，之后又将其推 广至三维问题。该方法在计算时，时间与空间步长无需满足c o u r a n t 稳定条件，因 此，与普通f d t d 相比，a d i f d t d 的时间步长可以增大至l o 倍左右，但数值色 散稍严重于普通f d t d 。该方法已广泛用于计算超大、超小尺寸目标的电磁特性。 另一方面，l i u 等人提出了一种时域伪谱算法。该方法利用了f f t 采样的无限精度 的特性，在空间上利用f f t 以及i f f t 完成对空间的微分。每个波长只需2 个网格。 2 第一章引言 从而使得计算量和存储空问大大下降。该方法同样在电大尺寸目标的散射问题中 得到了广泛的应用。 以上介绍了时域有限差分方法的基本发展过程。f d t d 方法结构简单，计算直 接，算法通用等优点，高阶f d t d 方法对数值色散也有良好的抑制作用，应用广 泛。研究高阶精度的f d t d 方法具有重要的学术意义和工程应用价值。 1 2 本论文的主要工作 本论文重点讨论了基于辛策略的高阶精度时域有限差分方法，同时做了如下 的创新性研究。 首先本论文推导了基于波动方程的高阶时域有限差分算法。该方法利用离散 奇异卷积方法对波动方程的2 阶偏微分进行差分近似，空域精度可达2 m 阶。同时 利用辛积分传播子方法构造4 阶时域精度。利用数值算例分析了波动方程形式的 高阶算法的精度特性，同时分析了该算法在计算复杂边界问题时的灵活性。以规 则波导的解析解作为参照，本论文的高阶算法比传统的时域有限差分方法的误差 低近2 0 0 d b 。在与传统时域有限差分方法达到相当的精度的情况下，高阶精度的算 法解决内存与计算量。同时计算了脊波导的t m 模式下的前几个截止波数，该结 果与已有文献的结果吻合。由此得出结论：基于波动方程的高阶算法有着良好的 精度特性，并且由于其简单的结构，对复杂边界条件的建模较为方便。 在已有的文献报道中，基于辛策略的高阶时域有限差分方法大多集中在处理 导波系统等问题中。对于处理电磁散射等开放边界类型的问题则较少有文献报道。 针对电磁散射问题，本论文推导了基于辛策略高阶时域有限差分方法的完全匹配 层的表示形式，利用修正的帕德近似公式对电导率和磁阻率进行近似，将高阶算 法中的p m l 反射误差进一步降低至已有文献报道的误差的1 。同时针对该算法 的特点修改了总场一散射场体系的连接边界条件。基于以上结论，利用基于辛策略 的高阶时域有限差分方法计算目标的电磁散射问题。计算了单个导体方柱、多个 导体柱的双站雷达散射截面。数值结果表明，在计算目标的散射特性时，高阶算 法与传统时域有限差分方法的计算结果吻合，但高阶算法使用的网格数低于普通 时域有限差分方法，从而有着更高的计算效率。同时，利用时域有限差分方法中 的周期边界条件，计算了一维金属带周期阵列结构的单个单元结构的散射特性， 计算该结构的单站和双站雷达散射截面。同时，进一步分析了一维金属带周期阵 列结构的吸收电磁波特性，我们计算了覆盖有该结构的金属导体柱的雷达散射截 电子科技火学硕士学位论文 面。结果表明，一维金属带阵列周期结构对电磁波有着一定的吸收效果。如何增 大前向散射，抑制后向散射，需要对该周期结构的结构参数和介质的特性参数做 进一步的调整和优化。 4 第二章时域有限差分方法 第二章时域有限差分方法 2 1时域有限差分方法的基本原理 时域有限差分方法是通过对m a x w e l l 方程中的微分直接进行差分近似得到的。 传统的f d t d 方法的精度为二阶。 在直角坐标系下，以其中一个方向为例，m a x w e l l 方程化为标量形式为： 其中盯和为介质的电导率和磁阻率。z 、为介质的磁导常数和介电常数。 若以图2 1 所示的y e e 元胞为基本单位对空间进行采样，以电场的工方向为例， 标量形式的m a x w e l l 方程可被离散为如下形式： e ( 嘎1 以尼) + z * 五1 + 秒1 ) 一矿* 三1 小五1 ，刁 a t f 1 s r j 日；+ ；( ，+ 三1 ，后+ 圭 一日；+ ：( ，+ j 1 ，后一芝1 ) l ya y ( 2 2 ) 按照以上的方式对所有的m a x w e l l 标量方程进行离散，按时间步进的方式进 行循环迭代运算，可得任一时刻电场或磁场的值。这种循环迭代运算是时域有限 差分方法的基础。 作为一种电磁场的数值计算方法，时域有限差分法具有一些非常突出的特点， 5 哦 q 怕 。 争等 。 邓 孚孚 出一0 、k-、一一尼 尼 一几 舭y一1 2 小一 0 卜可f可 兰扣磊砂化虿再不 一 + l l = 、 七 ， 1 2 + ， 卅 月x e 电子科技人学硕j ：学位论文 最主要包括以下几点： ( 1 ) 直接时域计算。它直接把含时间变量的m a x w e l l 旋度方程在y e e 氏网 格空间中转换为差分方程。在这种差分格式中每个网格点上的电场或磁场分量仅 与它相邻的磁场或电场分量及上一时间步该点的场值有关。它把各类问题都作为 初值问题来处理，使电磁场的时域特性被直接反映出来。 ( 2 ) 广泛的适用性。由于时域有限差分法的直接出发点是概括电磁场普遍规 律的m a x w e l l 方程，这就预示着这一方法应具有广泛的适用性。它可以模拟各种 复杂的电磁结构。媒质的非均匀性、各向异性、色散特性和非线性等均能很容易 地进行模拟。任何问题只要能正确对源和结构进行模拟，时域有限差分法就能给 出正确的解答。不管是散射、传输、透入或吸收的哪一种，也不论是瞬态问题还 是稳态问题。 ( 3 ) 计算程序的通用性。由于m a x w e l l 方程是时域有限差分法计算任何问 题的数学模型，因此它的基本差分方程对广泛的问题是不变的。此外，吸收边界 条件和连接条件对很多问题是可以通用的，而计算对象的模拟是通过给网格赋予 参数来实现的，因此一个基础的时域有限差分法计算程序对广泛的电磁场问题都 具有通用性。 ( 4 ) 简单直观，易于掌握。由于时域有限差分法直接从麦克斯韦方程出发， 不需要任何导出方程，这样就避免了使用更多的数学工具，使得它成为所有计算 电磁场的方法中最简单的一种。其次，由于它能直接在时域中模拟电磁波的传播 及其与物体作用的物理过程，因而它又是非常直观的一种方法。 6 第一二章时域有限差分方法 2 2时域有限差分方法的一些关键技术 2 2 1 稳定性及数值色散 在时域有限差分方法中，时间步长与空间步长不是相互独立的。为避免数值 解的不稳定性，必须满足一定的条件，即c o u r a n t 稳定性条件， c a t ( 2 3 ) 考虑波动方程的一般形式 磐a x + 磐o y + 孥o z 一丢c 尝a t = 。 c 2 4 ) 22222。 7 对于平面波 y ( x ，f ) = 厶e x p ( - 歹慨一纠) ) ( 2 5 ) 背+ 背+ 背铲射( 譬) 2 ( 譬) 2( 等) 2c 2 ( ! 警) 2 u。u ( 詈) 2 = + 蟛+ 碍 ( 2 7 ) 由公式( 2 6 ) 、公式( 2 7 ) 可知，当空间步长趋于0 时，数值色散可以减小 到任意程度。 2 2 2 截断边界条件 7 电子科技火学硕士学何论文 一般地，求解电磁场问题总是假定问题空间是无限大的，即所谓“开放 系 统。在采用f d t d 法计算时，所取研究空间越大，需要的存储量也就越大，计算 时间也就越长。由于计算机硬件条件的限制，实际上不可能模拟无限大的空间， 因此网格总要在某处被截断。不合理的截断对计算结果影响很大。因此希望通过 边界处理办法，使得空间的截断不影响或尽可能小地影响对实际问题的模拟。一 种行之有效的方法就是在截断处设置截断边界条件。所谓截断边界条件是这样的 一种算法：它不仅能保证边界场计算的必要精度，而且还能大大消除非物理因素 引起的入射到截断边界的波的反射，使得用有限的网格空间就能模拟电磁场在无 界空间中的传播。截断边界条件是f d t d 法的要素之一。截断边界条件的研究已 经有较长的时间了，目前这方面的进展较为迅速。最初t a y l o r 等人采用的是外推 法，效果欠佳且缺乏改进的一般理论。进入七、八十年代，出现了四大类截断边 界条件：基于s o m m e r f i e l d 辐射条件的b a y l i s s t u r k e l 截断边界条件；基于单向波 动方程的e n g q u i s t m a j d a 截断边界条件；利用插值技术的廖氏截断边界条件； m e i f a n g 超截断边界条件。这些截断边界条件通常在f d t d 仿真区域的外边界具 有0 5 5 的反射系数，在许多场合可视为无反射截断。进入九十年代，b e r e n g e r 和g e d n e y 分别提出的两类完全匹配层( p e r f e c t l ym a t c h e dl a y e r , p m l ) 截断边界条 件，可得到比上述截断边界条件低4 0 d b 的反射系数，使截断边界条件的研究又向 前迈了大步。 2 2 3 激励源设置 实际的电磁场问题总是包含有激励源，因此用f d t d 法分析问题时，应将被 研究的对象在真实源激励下这一完整条件在数值计算中尽可能近似地“复现出 来。在激励源的引入中，一方面需要考虑设置的源能尽量接近研究问题的真实激 励，另一方面需要考虑激励源的实现尽可能紧凑，以便尽量减少由此而来的计算 机内存占用和计算时间。激励源有多种，可以按照不同的型式和性质来分类。从 空间分布来看，有面源、线源、点源等。典型的面源有常见的平面波源。从频谱 特性来看，有工作与一个频率上的周期变化的连续波源，也有占较宽频谱的波源。 从源的时变特点来看，有随时间周期变化的时谐源，有对时间呈冲击函数形式的 波源，包括高斯脉冲、矩形脉冲、上升余弦脉冲等形式。 第三章高阶精度时域有限差分方法 第三章高阶精度时域有限差分方法 3 1 空域上的高阶精度算法 3 1 1 离散奇异卷积( d s o ) 方法 传统的f d t d 采用两点中心差分对时间和空间的一阶导数进行近似。由泰勒 展开式可知，这种近似所带来的误差是二阶无穷小量。如果增加插值的点数，误 差能够被控制在更小的范围内。邵振海等人提出了一种离散奇异卷积( d i s c r e t e s i n g u l a rc o n v o l u t i o n ，d s c ) 时域方法，其在计算电磁场问题时，引用d s c 方法离 散计算空间，构造高阶插值方程从而提高计算精度。下面简要介绍d s c 方法。 假设丁表示系统的特征函数，删为一个测试函数，则们的奇异卷积表示为 f o ) = 丁o ) 木厂o ) = l t ( , - x 旷( 沁k ( 3 1 ) 其离散的奇异卷积( d s c ) 形式为 v o ( , ) - - 瓦o 一沙k ) ( 3 - 2 ) 其中l i m 死g ) = r g ) ，规表示空间采样点。 若死为d s c6 核，由式( 3 2 ) 可得任一函数以及它的 阶导数为： 厂( ) g ) 以g 一以沙k ) ( 3 - 3 ) 其中万a ( 月, a ) & k - - x k ) = ( 丢r 以，g 一) ，疋a g 一) 即为巧核，m 是计算的带宽。 j 核可以是s h a n n o n6 核，亦可是l a g r a n g ej 核或p o i s s o n6 核。本文讨论拉格朗目 6 核( l a g r a n g ed e l t ak e m e l ，l k ) 的情况。与其它类型的6 核相比，l k 在相同的计 算带宽m 的情况下有更高的精度。 由于在实际的数值计算中，计算带宽不可能为无穷大，因此在必要时可以使 用以下的窗函数来限制公式( 3 3 ) 的带宽： 9 电子科技人学硕十学位论文 l k 的形式如下： 加窗后的l k 的形式为 陆唧( 一1 2 0 - 2 洲 4 ， k。阱k筒-mj，。网x-xj)y=km ，。g ) = 兀口 kv k “i j ( 3 5 ) 舡班恒k - m , i 。薏h 一譬 6 ， 对于波动方程，在式( 3 3 ) 中取n = 2 。在公式( 3 5 ) 与公式( 3 6 ) 中对x 求 二阶导数，取x e = o ，x f = x k + i a x ，当m = 6 时二者比较结果如图3 - 1 所示。 2 1 0 o 三 要1 - 2 3 图3 - 1m = 6 时未加窗的l k 与加窗的l k 的比较 由图3 - 1 可知加窗前后l k 变化很小，故窗函数不是必须的。 中n = 2 ，基于l k 的2 m 阶d s c 的误差项为 啪) = 销，蚤m 互m ，重g 1 ) c 为计算区间内的某个值。 3 1 2 离散奇异卷积方法的数值色散特性 对于公式( 3 3 ) ( 3 - 7 ) 在对空间偏微分的近似中，传统时域有限差分方法采用的是二阶中心差分， l o 第三章高阶精度时域有限差分方法 其精度为2 阶。 望l ：盈 也2 ( 3 - 8 ) 苏 ，缸 其中血为空间步长。假设f i - = j ( x i ) = e x p ( k x i ) ，带入式( 3 8 ) ，有 珐i n ( 三尬 ( 3 9 ) 左端为准确解波数= ，右端为波数的数值解，即。) = 2 s i n a k x ) 。两者的 逼近程度反映了差分方法的数值色散的大小。 图3 - 2 当带宽分别为1 ，2 ，3 时的数值色散 对于拉格朗日基的离散奇异卷积方法，当区单边计算带宽m = i 时，精度与传 统时域有限差分方法相当，对于m = i ，2 ，3 时，不同的数值色散关系如下所示。 冲s i n ( 去脸) ( 3 - 1 0 ) 缈( 缈) = 三s t n ( 三缈) 一五1s i n ( 三缈) ( 3 1 1 ) 电子科技火学硕士学位论文 缈( ) = 2 3 4 3 7 6 s i n ( 詈) 一。，3 。2 。8 4 s i n ( 三彩) + 0s i n ( 主缈) ( 3 - 1 2 ) 以准确解k a x 为横坐标，数值解为纵坐标，给出以上三种格式的误差比较， 如图3 2 所示。 由图3 2 可知，当离散奇异卷积方法的单边插值带宽m 增大时，波数的数值 解趋于准确解，即数值色散随着带宽的增大而减小。同时，随着空间步长的减小， 数值色散也会被抑制。 3 2 时域上的高阶算法 3 2 1 辛积分方法 首先考虑辛积分传播子。假设z 是一个向量，并且有 塑：c x( 3 1 3 ) d t 其中c 为一个矩阵。满足该方程的解x 在时间步长缸后可以由一个e x p ( a t c ) 来表 示，即 4 a 0 - - e x p ( m c ) x ( o ) ( 3 1 4 ) e x p ( a t c ) = i + a t c + 等c 2 + 等c 3 + 1 5 ) 其中缸o ) 为初值，且为单位矩阵。 如果c = a + b 且彳与b 不可交换，则e x p ( a t c ) 的值可以用辛积分传播子来逼近： e x p ( f c ) ：r ie x p p p 坍) e x p ( c p a t a ) + 。( 1 ) ( 3 16 ) 其中白和彩是表征传播子的实系数( 参见参考文献 1 6 】) 。n 为逼近的阶数，m 为 总的级数。 3 2 2m a x w e i l 方程的辛格式 考虑无源无耗各项同性的媒质中的m a x w e l l 方程的矩阵形式 1 2 第三章高阶精度时域有限差分方法 其中 k 一三d 肚e 甜肚 旦f ，h 1 ：形f ，日 o t l e jl e ，、 a u一 o z a o z a 砂 o a 砂 a o x a ，、 u 苏 ( 3 1 7 ) 日为磁场分量的向量，e 为电场分量的向量。，s 分别为介质的导磁率和介电常 数。 假设w = u + v u ：i0 3 x 3 l 0 ，。， f ，0 3 。30 3 x 3 、1 肚慨o ，。l s 由于阢矿不可交换，根据已有结论，可知 考虑到 e x p ( a f 形) ：ne x p p p a t v ) e x p ( e p a t u ) + 。( a t 州) ( 3 18 ) p = l f o x p ( a f u ) = 1 6 + 缸u l e x p ( f 矿) = 1 6 + a t v 因此e x p ( a t u ) ，e x p ( a t v ) 的值可以直接得到。 若 e 付日协 考虑二维t m 模式，二阶辛积分的形式为 1 3 ( 3 1 9 ) 、 r ，吆 电子科技大学硕十学位论文 彰+ 1 1 研一坐竽 印 日，n + ；：日，n + 一c t a t 竽 ( 3 - 2 0 ) 。 “o x 母竽p n + l - 雩 磁+ - ：风n + j 1 一_ _ c 2 a t _ a e - , 2 h 粤日，n + 三1 + 坐譬 ( 3 - 2 1 ) 。 na x 矽彰；+ 竽( 譬一等 占i僦卵j 其精度为2 阶。若要达到4 阶时间精度，则可以取阶数为5 。各阶数时的辛积 分系数如表3 - 1 所示。 表3 12 、4 阶精度的辛时域方法各级传播子系数 二阶精度 4 阶精度 二级传播子5 级传播子 c l o 50 1 7 3 9 9 6 8 9 1 4 6 5 4 1 矾 1 oo 6 2 3 3 7 9 3 2 4 513 2 2 c 2 o 50 1 2 0 3 8 5 0 4 1 2 1 4 3 0 杰 o oo 12 3 3 7 9 3 2 4 513 2 2 c 3 0 8 9 2 7 7 6 2 9 9 4 9 7 7 8 以 一 比 在辛积分中，各个传播子之间的关系为： 勺= + ，- p ，( o 掰+ 1 ) 彩2 如中，( o 印 ，1 ) 磊产0 对于一个2 级和5 级传播子的辛积分策略，其电磁场随时间步长推进的方式 如图3 3 所示。 1 4 第三章高阶精度时域有限著分方法 l 。r 。- _ _ l _ _ 。- - - - 。_ 1 。- - _ _ - 。- _ 。【 d h l t 2 h “h n ；1 n h ”l 夕l 呸巫萝恒巫砂| 呸剜l l 、 于 ， 屯矾i 乒亟聪乒 t 旷1 巧矗“t i ”1 1 5 h n + 狮矿粥 i r l 枷 小l 歹鞑剜恒亟刚恒剜羟垂画眶亟巫飚剜l 三一气， 喜专喜拿 r 专气， 越汹 商蒯崩乒 e ne 州肛e ”扔矿舢e 一( _ e 州勺 _ - _ _ - _ _ l _ - - - i _ - - _ _ _ - - - - _ - - - - _ i 上 。 取空间精度为4 阶，时间精度也为4 阶，t m 模式下的匠分量在n + 3 5 时刻的 值表示如下： e + 3 巧( f ，) = e + 2 艏( f ，) + 警曙渺5 ( 1 ，小矿5 ( 气1 ，切 一上2 4 e 伽。y 俗( 一) 掣巧( 叫3 ) ) 协2 2 ) 3 3 波动方程的高阶时域有限差分格式 各向同性的无源媒质中电场波动方程为 v 2 脚占警 其中小分别为介质的导磁率和介电常数，以及 ( 3 2 3 ) 电子科技人学硕十学位论文 v 2 = 筹+ 熹+ 轰a x za v z貌2 如果设豆：望为中间变量，则有 西 貉睁) f ，瓯， 矿：iv 2 ， l 石厶 ) 0 3 。3 - - u + z u ：f ，q ， l q 。， f ，q 。， 矿：iv 2 ， i 一3 l s 由于队y 不可交换，应用s i p 方法，得到 ( 3 2 4 ) ( 3 2 5 ) ( 3 2 6 ) ( 3 2 7 ) ( 3 2 8 ) e x p ( a t i v ) ：ne x p 0 ，a t v ) e x p ( c p a t u ) + 。( a t 肘- ) ( 3 2 9 ) 7 = 1 同样有 e x p ( a t u ) = 6 + a t u ( 3 - 3 0 ) l e x p ( a t v ) = 1 6 + a t v 对于四阶精度的s i p 方法，m = 5 。c l ，c 5 ，d l ，d 5 的值参见表3 1 。 根据以上的分析，波动方程的高阶算法实际需要存储2 个分量，即公式 ( 3 2 5 ) 中的e 和易。而传统的y e e 方法在计算t e t m 问题时，至少需要3 个分 量。因此本论文的方法可节约将近1 3 的内存。同时，由于该方法避免使用y e e 网格，因而减少了计算的复杂性。 1 6 一、 、_、 姆蝎 厶眇 凸伪 第三章高阶精度时