（电磁场与微波技术专业论文）雷电电磁脉冲的时域有限差分法研究.pdf

摘要 雷电是自然界的放电现象。雷电发生时，在极短的放电时间里能产生极高的放电电压， 其间释放出的巨大电能可以在几十微秒内把雷电通道加热至1 0 0 0 0 k 以上，同时产生的强 烈的电磁辐射一雷电电磁脉冲( l e m p ) 会对周围的电力系统、计算机和各种电子设备造成干 扰甚至瘫痪。在当今这个信息社会里，雷电被称为是电子时代的一大公害，国内外越来越 多的专家在着手进行l e m p 方面的理论研究和仿真计算。 本文的主要研究内容包括：介绍了雷电的基本物理过程：作为计算工具介绍了时域有 限差分法( f d t d ) 的基本原理，研究了m u r 吸收边界条件；讨论了常用的几种雷电通道 底部电流函数并画出其波形图；介绍了几种雷电回击模型，重点研究了最常用的工程模型， 并对各种模型的应用条件及优缺点进行了讨论；以雷电通道底部电流为已知条件，采用各 种回击模型，使用f d t d 方法计算了雷电通道周围的电磁场。首先以简单高斯脉冲函数作 为通道基部电流函数，把大地当作理想导体处理，采用传输线模型( t l ) 计算了雷电电磁 脉冲，并和文献 1 1 的结果做了对比；接着采用双指数和h e i d l e r 函数的叠加函数作为通 道基部电流函数，把大地当作有耗介质处理，分别采用传输线模型( t l ) ，线性衰减的传 输线模型( m t l l ) ，指数衰减的传输线模型( m t l e ) ，b g 模型，t c s 模型计算了通道周 围的电磁脉冲，并给出电磁脉冲的时域波形图，计算结果与文献 2 进行了比较，得出一 致的结果；最后把雷电通道当作理想导体处理，受通道底部一实际的电流源激励，在f d t d 计算中采用细导线的特殊处理方法，计算出了雷电通道中不同高度处的电流分布及周围的 电磁场，并把这种细导线模型和t l 模型的结果进行了比较分析，得出了和文献 1 1 一致 的结论。 关键词：雷电电磁脉冲，时域有限差分法，回击模型，通道基部电流，细导线模 型 i i a b s t r a c t l i g h t n i n gd i s c h a r g ei s o n eo fn a t u r a lp h e n o m e n o n s t h e r ei s v e r y h i g hv o l t a g ew h e n l i g h t n i n gd i s c h a r g e ，t h ee l e c t r i c a le n e r g ed e p o s i t e di n t ot h ec h a n n e li nt h ef o r mo fh e a ta n dt h e l i g h t n i n gc h a n n e li sh e a t e dt ot e m p e r a t u r ea b o v e10 0 0 0 ko rs ow i t h i naf e wt e n so fu s t h e e l e c t r o m a g n e t i cf i e l dg e n e r a t e db yl i g h t n i n gi sr e s p o n s i b l ef o rd a m a g e so rd i s t u r b a n c e so f d i f f e r e n tk i n do fe q u i p m e n ts u c ha sp o w e rt r a n s f o r m e r s ，c o m p u t e r sa n di n g e n e r a le l e c t r o n i c d e v i c e s w i t ht h ec o m i n go fi n f o r m a t i o ne r a , t h ea c c i d e n to fd a m a g et oe q u i p m e n to rs y s t e m c a u s e db yl i g h t n i n ge l e c t r o m a g n e t i cp u l s e ( l e m p ) i si n c r e a s e dh i g h l ya n dm a n yl i g h t n i n g r e s e a r c h e r so fa l lo v e rt h ew o r l d p a ym o r ea t t e n t i o nt ol e m p t h ec o n t e n to ft h i st h e s i si so r g a n i z e da sf o l l o w i n g ：i nc h a p t e ro n e ，t h el i g h t n i n gp h y s i c a l p r o c e s s e sa r ei n t r o d u c e d i nc h a p t e rt w o ，t h eb a s i cp r i n c i p l eo ff 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 ) i sp r e s e n t e da n dt h em u ra b s o r b i n gb o u n d a r yc o n d i t i o ni ss t u d i e d i nc h a p t e rt h r e e ， f i r s t l yaf e wo fm e t h o d st om o d e lt h el i g h t n i n gr e t u r ns 缸o k ea r ei n t r o d u c e da n dt h e ns o m e e n g i n e e r i n gm o d e l sa r ei n v e s t i g a t e d ；s e c o n d l ys e v e r a lc u r r e n tf u n c t i o n so fr e t u r ns t r o k ec h a n n e l b a s ea r ei n t r o d u c e da n dg r a p h sa r em a d ef o rc o m p a r i n gw i t ht h er e p r e s e n t a t i v er e t u r ns t r o k e c u r r e n tm e a s u r e dp r a c t i c a l l y t h i r d l yl e m pr a d i a t i o ni sc o m p u t e d 、i t hf d t dm e t h o d i nt h e f i r s tc a l c u l a t e de x a m p l e ，b ya s s u m i n gt h el i g h t n i n gc h a n n e lt ob es t r a i g h ta n dv e r t i c a la n dh a sn o b r a n c h e sa b o v ep e r f e c t l yc o n d u c t i n gg r o u n d ，w eu s e dt h es i m p l eg a u s sp u l s ea st h ec u r r e n t f u n c t i o no ft h er e t u r ns t r o k e ，t h et r a n s m i s s i o n - l i n e ( t l ) m o d e li su s e dt oc o m p u t e dl e m p t h e n w eu s eh e i d l e rf u n c t i o nt os i m u l a t ec h a n n e lb a s ec u r r e n t ，t h eg r o u n dh a saf i n i t ec o n d u c t i v i t y ， t h el e m pr a d i a t i o n i sa l s oc o m p u t e d i nt h i sc o m p u t a t i o nw eu t i l i z e dt h et r a n s m i s s i o nl i n e m o d e l ( t l ) ，m o d i f i e dt r a n s m i s s i o n l i n em o d e lw i t hl i n e a rc u r r e n td e c a yw i t hh e i g h t ( m t l l ) ， m o d i f i e dt r a n s m i s s i o nl i n em o d e lw i t h e x p o n e n t i m c u r r e n t d e c a y w i t h h e i g h t m e t h o d ( m t l e ) ，b r u c e g o l d em o d e l ( b g ) a n dt c sm o d e l ，a n di l l u s t r a t e dt h ed i s t r i b u t i o no f e l e c t r i cf i e l d a tl a s t ，a s s u m i n gt h el i g h t n i n gc h a n n e lt ob eap e r f e c t l yc o n c u c t i n gv e r t i c a lw i r e e x c i t e db yal u m p e dc u r r e n ts o u r c ea ti t sb a s e ，w h o s ee l e c t r o m a g n e t i cf i e l ds t r u c t u r ei sa n a l y z e d u s i n gf d t d a l lc a l c u l a t i o ni sc o m p a r e d a n dt e s t e dw i t ht h er e s u l to f p u b l i s h e dl i t e r a t u r e k e y w o r d s ：l i g h t n i n ge l e c t r o m a g n e t i cp u l s e ( l e m p ) ，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 ) ，r e t u r n s t r o k em o d e l ，c h a r m e lb a s ec u r r e n t ，t h i nw i r em o d e l i i i 南京邮电大学学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究 工作及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的 地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包 含为获得南京邮电大学或其它教育机构的学位或证书而使用过的材 料。与我一同工作的同志对本研究所做的任何贡献均已在论文中做了 明确的说明并表示了谢意。 研究生签名：牲日期：j 垒尘辱二必 南京邮电大学学位论文使用授权声明 南京邮电大学、中国科学技术信息研究所、国家图书馆有权保留 本人所送交学位论文的复印件和电子文档，可以采用影印、缩印或其 它复制手段保存论文。本人电子文档的内容和纸质论文的内容相一 致。除在保密期内的保密论文外，允许论文被查阅和借阅，可以公布 ( 包括刊登) 论文的全部或部分内容。论文的公布( 包括刊登) 授权 南京邮电大学研究生部办理。 研究生虢肛新虢嵇年魄也阻哲 南京邮电大学硕士学位论文 第一章绪论 1 1 研究背景 第一章绪论 人类对于闪电的最早研究是在1 8 世纪中叶由本杰明富兰克林( b e n j a m i nf r a n k l i n ) 所完成的。两百多年前，富兰克林在雷雨中放风筝，由他冒着生命危险而得出的实验成果， 证实闪电是带电的。富兰克林在后来的实验中发现，若将一跟长铁杆与地面接触，铁杆在 雷击之后会静静地将雷电导引至地底下去，顶多只出现一些火花。这种现象与物体遭雷击 相比，实在温和太多，于是乎有避雷针的发明，使得人们可借以避免雷击。虽然如此，人 类仍然年年因闪电而蒙受巨大损失。在富兰克林之后，闪电的研究一直没什么进展，直到 1 9 世纪末，当照相术( p h o t o g r a p h y ) 和光谱学( s p e c t r o s c o p y ) 成为研究闪电的利器后，人 们对闪电才有更进一步的了解。现代对闪电的研究开始于2 0 世纪，以c t r 维尔森在该 领域所做的研究工作为标志。维尔森进行了远程闪电电场测量与分析，因而成为第一个能 够推测雷雨云电荷结构，以及闪电时所释放的电量的人。雷电是联合国确定的全球最严重 的十种自然灾害之一。全球每年因雷击造成的人员伤亡、财产损失不计其数，导致火灾、 爆炸、信息系统瘫痪等事故频繁发生，从卫星、通信、导航、计算机网络乃至每个家庭的 用电设备都会受到雷电灾害的严重威胁。雷电造成的灾害可分为直接效应和间接效应。所 谓直接效应是指在雷击点处由雷击直接造成的破坏，除直接雷击外，还有雷电的静电感应 作用，闪电放电时的电磁感应作用；闪电放电时产生的强烈电磁脉冲；雷电反击以及雷电 过电压波可能沿各种架空线、无线电天线、天线馈线、电缆外皮和金属管线等传入仪器设 备，酿成祸患。随着现代信息网络技术和微电子技术的高度发展以及在通信、电力、航空 航天等领域的广泛应用，使雷害对象也发生了转移，从对建筑物本身的损害更多地转移到 对电力系统、电器、电子设备的损害上来。随之防雷对象也由强电转移到弱电。雷电产生 的电磁感应已成为当今对人类的主要危害。研究雷电放电时产生的雷电电磁脉冲 ( l i g h t n i n ge 1 e c t r o m a g n e t i cp u l s e ，l e m p ) 已经成为当今雷电防护领域的一个重点。 随着高科技的发展，雷电灾害会显得越来越严重，所以对l e m p 做深入的研究就显得非常 必要和迫切。 1 2 雷电的物理过程 闪电通常在雷雨云情况下出现，但在雨层云、雪暴、尘暴和火山爆发时也会偶尔出现。 闪电按发生的部位可分为四种型态：第一种是云层内部的放电现象( i n t r a c l o u d 1 南京邮电大学硕士学位论文第章绪论 d i s c h a r g e ) ，占所有闪电的绝大部分；第二种为云对周围空气的放电现象( a i r d i s c h a r g e ) ，通常发生在云顶：第三种是云与云间的闪电：最后一种是云对地面的闪电 ( c l o u d t o g r o u n d1 i g h t n i n g ) ，这是最广为研究的类型，因为它们对人们的生命财产有 极大的威胁性。前三种闪电统称为云闪。本文对闪电的研究主要是针对第四种，与人类生 活息息相关而且也是被研究最多的类型云对地的闪电，简称为地闪。地闪根据所中和 的云中电荷极性又可分为正地闪和负地闪，云层带正电荷对地放电称为正闪电，而云层带 负电荷对地放电称为负闪电。一般情况负地闪约占总地闪的9 0 。 闪电，人眼看上去似乎是一次瞬间闪光，但通过高速摄影揭示，它往往是由同一条通 道、彼此间隔约百分之几秒的多次相继放电组成( 图1 1 ) 。“一个闪光”( af l a s h ) ， 这是一整个放电过程( t o t a l d i s c h a r g e ) ，时间约可持续0 2 秒，一次闪电是由多个闪光所 组成。一个闪光又可由三至四个称为“一击”( as t r o k e ) 的放电单元所组成，每一击持 续的畴间约为几十个毫秒( m i l l i s e c o n d ) 。而每一击是由一个光度较弱的先遣放电过程 ( p r e d i s c h a r g e ) ，称为“先导过程”( 1 e a d e rp r o c e s s ) 与伴随而来的回击( r e t u r n s t r o k e ) 所构成。先导过程是由云层传递至地面；先导过程之后是一个快速的、高亮度由地面传 播至云层的回击( r e t u r ns t r o k e ) 。每次闪光的第一击，其先导过程称为“阶梯先 导 ( s t e p p e d l e a d e r ) 。阶梯先导开始于雷云下半部之负电区与底部的零星正电区之间局 蛳雠 r 一一2 0 心e c 1 r ，o 图1 1 闪电的过程，时间往右增加。【摘自m a u m a n “l i g h t n i n g ”( 1 9 6 9 ) 】 2 til要；：fll上 南京邮电大学硕士学位论文第一章绪论 部的电崩溃( e l e c t r i c a lb r e a k d o w n ) ，这使得原先存在于冰粒或小水滴上的电荷更加活泼 ( m o b i l e ) ，最后云底的负电荷产生了足够的电场，可以开辟一条通路到达地面，这个通路 称为阶梯先导。每个阶梯( s t e p ) 约长5 0 公尺，阶梯与阶梯之间，间隔约5 0 微秒 ( m i c r o s e c o n d ) ，接着下一个阶梯产生。阶梯先导一般的平均速率为1 5 1 0 5 公尺秒， 约为光速的两千分之一。若云底的高度是3 公里，那么阶梯先导从云到地面将花2 0 毫秒 的时间，其平均电流约为数百安培，这路径的半径约在1 到1 0 公尺，将传递约5 库仑的 电量至地面。回击( r e t u r ns t r o k e ) 是跟在阶梯先导之后，当阶梯先导之通路接近地面就 像放了一根导线，强大的电流以极快的速度由地面流至云层。这一个过程，称为回击， 约需7 0 微秒的时间。回击的平均速率约为光速的三分之一至十分之一，典型的回击电流 强度约为一至两万安培。回击的发光强度比先导强得多，肉眼所见的闪光即为回击。 回 击的亮度相当耀眼，电流很快达到尖峰值( 所需时间小于一毫秒) 并持续几毫秒，再落到峰 值的一半，并再持续2 0 至6 0 微秒，然后电流减弱到几百安培，并持续几个毫秒。回击 的通路( c h a n n e l ) 温度可高达约3 0 0 0 0 。k ，高温使得通路中的气体膨胀，于是使得通路急 速膨胀，其膨胀速度超过音速，于是产生了音爆，这就是我们听到的雷声( t h u n d e r ) 。经 过约5 至1 0 微秒，通路中的空气与周围空气达成平衡，此时通路的直径缩到约几公分大 小。如果雷电的电流就此停止流动，那么闪电之闪光也就结束，但若雷云仍有放电的能力， 将继续有第二击或第三击。第一击之回击结束后，当雷云又拥有足够的电场时，第二击 即开始。第二击开始与上次回击结束间隔约几十毫秒。第二击的先导过程称为迅速先导 ( d a r tl e a d e r ) 。当回击与迅速先导之间隔小于1 0 0 毫秒，迅速先导将循着上一次回击的 路径。当间隔时间大于1 0 0 毫秒时，迅速先导会另辟路径，但路径上的某些点仍会变换 到原有的回击路径上；当间隔时间大于几百毫秒，那就是另一个闪光的阶梯先导了。如 同阶梯先导，迅速先导的阶梯长度也大约是5 0 公尺。传播速率比阶梯先导快，速率约为 2 10 6 公尺秒。 而迅速先导传递至地面的电量也比阶梯先导少。如同阶梯先导，在迅 速先导之后，接着第二次子回击。当雷云仍有足够电场可以放电，将会有下一个迅速先导， 若雷云已没有能量进行下一次放电，此次的放电过程就结束，这整个放电过程就称为一个 闪光( af l a s h ) 。一个闪光通常包含三至四击，而一次闪电将包含多次的闪光。 云闪的危害主要表现在航空、航天方面。而地闪，由于其产生的回击电流幅度很大， 一般为几十千安，甚至可高达上百千安，其产生的l e m p 危害更为显著，是目前l e m p 研究的主要方向。 南京邮电大学硕士学位论文第一章绪论 1 3 雷电电磁脉冲建模方法 在l e m p 理论计算方面，雷电回击模型的建立，使得人们可以通过数学计算对地闪 的辐射电磁场有很好的认识。1 9 4 1 年b r u c e 和g o l d e 1 】首次提出了雷电回击模型。此后， 国内外雷电研究者分别从不同的角度提出了各种雷电回击模型。1 9 9 8 年r a k o v 和u m a n 【2 1 对各种模型做了归纳，主要可以分为气体动力学模型、电磁模型、分布电路模型、工程模 型。我们将在后面的第四章详细介绍这些不同的雷电回击模型。 雷电模型建立起来，人们就可用它来研究计算电磁脉冲，雷电电磁脉冲的计算包括解 析方法和数值计算方法。 解析方法一般假设雷电回击通道无分支且垂直于地面，通道周围为无穷空间，通道电 流分布已知。这时通道中某- d , 段电流元i d h 可以看作是一个垂直电偶极子，我们可以求 得该偶极子在有耗上半空间的电矢位烈；。： 式中： 地= 警 兰掣+ 甓掣一2 r 厶c 桫舭删焘刊 民，r ，分别是观察点距电偶极子及其镜像偶极子的距离。 0 1 ，s l = c o e ，分别为大地电导率和电容率。 砖= 厢，i ：0 ，1 。 恐：扩砑；r ，：扩丽 醢= 一；w 。峙。+ ；w 硖= 一j w p 。b 。+ j w s 3 4 ( 1 1 ) 南京邮电大学硕士学位论文第一章绪论 由电偶极子产生的场可按下式计算： 皿= 而1 石面d 2 a ：oj w p o s od r d z ( 1 2 a ) 扭：= j w , u l o s o ( - 再d 2 z + 七； 彳：。 c1 2 b ) 吗= 一瓦1 了d a z o ( 1 2 c ) 式( 1 1 ) 中的第三项既是s o m m e r f e l d 积分，要求得s o m m e r f e l d 积分的解需要进行 近似处理，围绕着s o m m e r f e l d 积分就产生了各种近似方法，如理想近似法1 3 1 最速下降 法【4 】【5 1 、n o r t o n 6 】【7 】法等。 常见的用来计算雷电回击的数值方法主要有有限元法( f i n i t ee l e m e n tm e t h o d ，f e m ) 、 毪= l - t 自由空间嘞伊觎颤潮确 i ；j r - l 。v 。 、 旭。 驴 t r r i ： h 。 歹彳一_ 基磊。! 薹 j l _ i i i 州 i i l l 盂： - 图1 2 雷电回击通道模型 5 南京邮电大学硕士学位论文第一章绪论 矩量法( m e t h o do fm o m e n t s ，m o m ) 、时域有限差分法( 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 ) 。s h o o r y 和m o i n i 等人【8 l 把雷电的回击通道看成一个由底部电流源激励的、 垂直于地面的天线，使用m o m 方法在频域内求解电场积分方程，作为结果的通道电流分 布可以用来计算距离通道任意距离的电磁场分布。文献 9 】使用有限元法求解描述雷电回 击通道上的瞬时电压和电流分布的传输线方程，得出了距离回击通道远处的电压一时间和 电流一时间的波形。文献 1 0 1 提出了解析方法和f d t d 方法混合的方法计算辐射电磁场， 首先使用镜像法计算图1 2 所示偶极子的磁场分布，通过对式( 1 2 c ) 积分，求得f d t d 计算空间中网格点的磁场分布，然后使用f d t d 方法进一步求出回击通道周围空间的电 场分布，这个方法的提出由于对于回击通道近场的电场分布，使用不同的模型计算，不同 文献得到结果有很大的差别，然而对于磁场的结果却有一致的结论。文献 1 0 】的方法解决 了这个问题。文献 1 1 贝j j 完全使用了f d t d 的方法计算了回击通道上电流分布和近场的电 磁场分布。 1 4 本文的主要研究内容 本文主要研究地闪雷电在回击过程中产生的雷电电磁脉冲( l e m p ) 的数值计算，重 点是基于f d t d 方法对各种回击模型进行比较研究。 雷电是一个瞬间发生的过程，对于雷电电磁脉冲的计算，传统的解析方法是围绕 s o m m e r f e l d 积分进行近似处理，各种近似方法包括：理想近似法【3 1 、最速下降法【4 1 【5 1 、 n o r t o n l 6 1 1 7 1 法等。数值计算方法一般有有限元法【9 】、矩量法【8 】、时域有限差分法【1 0 】【1 1 1 等， 本文使用时域有限差分法( f d t d ) 基于各种回击模型计算了雷电回击通道的辐射电磁场。 主要计算了距离回击通道5 0 米、1 0 0 米处的电场分布，并和相关文献的结果进行了比较。 数值方法计算回击电磁场脉冲可以避免各种数学近似，不需要进行繁琐的积分、微分计算， 但是由于时域有限差分法受计算时间和计算空间的影响，所以在计算通道附近的电磁场分 布时比较方便。对于距通道较远的电磁场分布，比如距通道几千米远，f d t d 方法的使用 就受到一定的限制。本文在计算回击辐射电磁场之前还对f d t d 方法进行了全面的介绍， 并对几个经典的算例进行了计算。 6 南京邮电大学硕士学位论文第二章f d t d 算法理论 第二章f d t d 算法理论 2 1 时域有限差分法( f d t d ) 理论 k s y e e 1 2 1 于1 9 6 6 年提出了用空间网格和时间步的耦合偏微分方程组代替空间和 时间连续的耦合偏微分方程组( m a x w e l l 旋度方程组) ，奠定了用时域有限差分法 ( f d t d ) 计算电磁波问题的基础。g m u r 1 3 】于1 9 8 1 年提出了网格区域截断边界的 吸收条件( a b c ) ，进一步扩展了f d t d 法分析电磁波问题的范围。近2 0 年来，经过 一批科学家的不断改进，f d t d 方法逐步走向成熟，解决了一系列的问题： 吸收边界条件的应用和不断完善。由于用f d t d 求解问题一般假定问题空间是无限 大的，即所谓是“开放 系统。所以边界条件的设置很重要。在y e e 最初方法中使 用的是硬边界，即把边界设置为理想导体，这样就不可避免的在边界处引起电磁波 的反射。后来c t m u r 给出了一阶、二阶近似吸收边界条件，才使得边界问题得以不 断完善。 总场区和散射场区的划分。f d t d 在应用于散射问题的时候，有必要把入射波和散 射波加以分离，发展后的时域有限差分法利用连接边界条件把计算网格空间划分为 内部的总场区和外部的散射场区，吸收边界条件用在散射场区。 亚网格技术。传统的时域有限差分法中对物体的模拟的最小尺度是一个网格，为了 模拟细微的结构只能把网格划细，这样会大大增加对存储空间和计算时间的要求。 所谓亚网格技术就是模拟小于一个网格尺度的结构的方法。这样就扩大了时域有限 差分法应用范围。 广义正交曲线坐标系中的差分格式和非正交变形网格。这种方法把网格的形式多样 化，便于更精确的模拟各种形状的物体。 适用于色散介质的差分格式。在用时域有限差分法计算脉冲电磁场与色散介质的相 互作用时，由于脉冲含有宽的频率成分，差分格式中必须考虑介质的色散性质，传 统的差分格式己不再实用。近来用迭代技术解决了卷积计算的问题，从而获得了适 用于色散介质的差分格式。 超吸收边界条件和色散边界条件。传统的边界条件只考虑边界上的切向电场或切向 磁场，超吸收边界条件则既考虑电场需要满足的条件，又通过磁场的计算进行修正， 从而进一步提高了相应吸收边界条件的效果。色散吸收边界条件是为了适应不同频 率时传播速度不同的情况。 7 南京邮电大学硕士学位论文 第二章f d t d 算法理论 2 1 1m a x w e l l 方程 f d t d 算法就是用有限差分式替代m a x w e l l 方程对时间和空间的微分而得到的，所以 首先我们要了解以下包含有不同形式媒质参数的m a w e l l 方程。 对于任意媒质，m a x w e l l 方程的微分形式是： v h ：望+ ， 西 ( 2 1 ) v 小一詈 ( 2 2 ) b = 日 ( 2 3 ) d = 占e( 2 4 ) v b = 0( 2 5 ) v d = p ( 2 6 ) 以上式中e 是电场强度( v m ) ；d 是电位移( c 朋2 ) ；h 是磁场强度0 肌) ；b 是磁 感应强度( t ) ；p 是电荷密度( c m 2 ) ；，是电流密度( a m 2 ) 。 同时了解e 和d 及h 和b 间的关系是很重要的。在真空或自由空间中，它们间存在 简单的比例关系，即： b 。心日 ( 2 7 ) d = c o e 对于各向异性媒质，上述关系变为下式： b 2 日 ( 2 8 ) d = g e 式( 2 7 ) 中, g o ，毛分别为自由空间中的磁导率( 日所) 和介电常数( f m ) ，式( 2 8 ) 中，s 分别为媒质中磁导率( 日历) 和介电常数( f 肌) ，, g = l o g ，g = 氏占，。g r , ，分别为相 对磁导率和相对介电常数。 2 1 2f d t d 基本方程 先将m a x w e l l 旋度方程在直角坐标系中展开为六个标量场分量，再将问题空间沿三 个坐标轴向分成很多网格单元，每个网格单元长度作为空间单元，相应得出时间单元。用 南京邮电大学硕士学位论文第二章f d t d 算法理论 有限差分格式表示关于场分量对空间和时间变量的微分，即可得到f d t d 基本方程。 我们假定研究的空间是无源的，并且媒质参数s ，u ，仃，s 不随时间而变化，在坐标系 ( x ，y ，z ) 中，m a x w e l l 方程将化为以下六个标量方程： ( 2 9 ) 式( 2 9 ) 中占，仃，s 分别为媒质磁导率( 日m ) 、介电常数( f 研) 、电导率 聊) 和磁电 阻率心i n ) 。 上述六个方程构成了完整三维问题的情形。在实际中，往往会遇到研究的模型沿一个 轴向或两个轴向不变化的问题，对于前者是一个二维问题，后者就构成了一个一维问题。 对于二维或一维问题只要把上述方程加以简化就可以。 图2 1基本空间单元上场分量图 如图2 1 所示，用一长方体将电磁问题的求解域包含在内，并溉只z 三个方向将该长 9 e 髟 t a一钟a一所a一研 8 占 s + + + 盯 盯 仃 = = = 哆一瑟啦一缸啦一砂 一 一 一 旦砂盟瑟堡瑟 r r = 以 耳 见 a一钟矿一彭a一国 心 知 心 + + + o 0 0 i i = = 诅一钞甄一瑟哆一(毽 一 一 一 堡瑟堡叙堡钞 南京邮电大学硕士学位论文第二章f d t d 算法理论 方体用网格离散，网格步长分别为位，少，止，网格节点的标号分别以f ，猿示。这 时，第( f ，助个节点的坐标可表示为： ( 薯，咒，刁) = ( i a x ，j a y ，k a z ) ( 2 1 0 ) 若将时间轴也以时间步长进行离散，则在第外节点上第n 个时间步的任一场量值可以表示 为： f ”( f ，尼) = f ( i a x ，f 缈，尬，n a t ) ( 2 1 1 ) 这里f 表示任一场量，、，和刀为整数。为了实现关于空间坐标与时间变量的差分近 似，并考虑到电磁场在空间相互正交和铰链的关系，y e e 提出了如图2 1 所示的差分网格 单元。考虑到y e e 网格中六个场分量的相对位置和表示式( 2 9 ) ，将这六个场分量所满足的 一阶耦合偏微分方程组( 2 9 ) 中关于空间和时间变量的偏导数用中心差商近似，则可得 如下时域有限差分方程组： 爿：+ ( f ，歹+ 吉，后+ ) = 月：一i ( i ，+ 专，后+ 吉) 。 出 i 髟( f ，+ ，忌+ 1 ) 一髟( f ，_ ，+ 吉，七) e ( f ，j + 1 ，后+ ) 一霹( f ，七+ ) l ( f ，+ 圭，k + ) i 心 缈 l ( 2 1 2 ) 月哆+ 5 ( f + 吉，七+ 士) = 爿哆一5 ( f + 专，后+ 号) + 一 ( f + 吉，j ，七+ 圭) ie ( f + 1 ，j j + 吉) 一e ( f ，尼+ ) e ( f + ，后+ 1 ) 二e ( f + ，后) l l a xa z l ( 2 1 3 ) 彤+ 5 ( f + 吉，歹+ 专，尼) = 彰一5 ( f + 壬，+ 号，尼) 。 出 lp ( f + 吉，+ 1 ，后) 一e ( f + 圭，后) e ( f + 1 ，+ 专，七) 一髟( f ，+ ，后) i ( f + ，_ ，+ 吉，k ) i a y a x i i o ( 2 1 4 ) 南京邮电人学硕士学位论文 第二章f d t d 算法理论 f “( ，+ ，_ ，后) = 主三g l ：三 i 囊詈i 黼fo + 吉，七) + 乏；石j _ f ：i 乏了2 ：a ：t i 了：丽 l 笙：垫i ：i ：查! 二笙：垫i ：二i ：盟一笙：垫i ：生i ! 二垡：垫i ：苎二望l i缈 垃 l ( 2 1 5 ) 髟“(l+，七)=主三gi手圭毒妻i糕髟(f，+，七)+五：石：f_ij弓j2ia；ti：丽 ( 2 1 6 ) e“(，七+)=主要罢乡l圭i端f(，七十吉)+三：石_zl_三i2；a：t五z了面 i 堡：堑i ：z ：生i ! 二竺：垫二i ：生二望一笙：垫i ：生i 2 二丝：垫二i ：生兰l i 缸 缈 i ( 2 1 7 ) 由于关于空间和时间变量的偏导数都采用了中心差商近似，以上时域差分方程的收敛 阶数为o ( a u 2 , a t 2 ) ，a u = m a x ( a x ，a y ，a z ) 。 从方程可以看出，空间网格节点上某一时间步的电场值取决于该点在上一时间步的电 场值和与该电场正交平面上相邻节点处在上半时间步上的磁场值，以及媒质的电参数o 和 ；空间网格节点上某一时间步的磁场值取决于该点在上一时间步的磁场值和与该磁场正 交平面上相邻节点处在上半时间步上的电场值，以及媒质的磁参数u 。 2 1 3 边界条件 将时域有限差分法应用于开放区域上的电磁问题时，需要将开放域截断为有限区域。 这时，在截断边界上必须引入截断边界条件( t e r m i n a t e db o u n d a r yc o n d i t i o n ) 或吸收边界 条件( a b s o r b i n gb o u n d a r yc o n d i t i o n ) ，以模拟被截去的外部空间的影响。一个好的吸收边 界条件应在截断边界非常靠近物理结构不均匀区域时仍然能获得正确的满足精度要求的 解。截断边界非常靠近物理结构不均匀区域意味着需要进行网格剖分的求解区域减小或者 说网格节点数的减小，从而有效地减少所需计算机储存空间和计算时间。因此，吸收边界 条件是时域有限差分法的一个重要方面。国际上近年来对吸收边界条件的卓有成效的研究 南京邮电大学硕士学位论文第二章 f d t d 算法理论 有力地推动了时域有限差分法的发展。比较著名、应用广泛的边界条件有：基于单向波动 方程的e n g q u i s t - - m a j d a 吸收边界条件n 钉：使用插值技术的廖氏吸收边界条件n 射；梅一方 超吸收边界条件n 钔；b e r e b g e r 完全匹配层n 铂等。在这些边界条件中，m u r 给出的单行波 方程的差分格式比较简单适合于柱坐标f d t d 网格，本文仅仅介绍m u r 吸收边界条件。 假设只考虑自由空间中的正立方体网格，用c 表示电磁波在真空中的光速，骧示电场 和磁场六个分量中的任何一个。于是啪满足以下标量方程 1 ( a ：+ a ；+ a ；一击a ；) f = 0 ( 2 1 8 ) f 考虑网格截断面萨0 处的吸收边界条件，其余网格处于x 0 。将上式作算子分解，有： a 。一a ，】 a ，+ a ， 上式又可以分解为： ( a ，- v o ，) f = 0 ( a 。+ 旧，) f = 0 其中，= 方程的解为：f ( x ，y ，z ，f ) = 甲o + v x ) 】f = 0 ( 2 2 1 a ) f ( x ，y ，z ，) = 甲p w ) ( 2 2 1b ) 其中第一个方程描述的是沿吱方向传播的行波，而第二个方程描述的是沿h 方向传播的行 波。 。 理想的吸收边界条件意味着在截断边界上没有反射，换句话说在x = d 处应只有沿叫方向 传播的波而没有沿h 方向传播的波。由于方程( 2 2 0 a ) 的解( 2 2 1 a ) 就是一个沿曩方向传播的 行波，因此将方程( 2 2 0 a ) 应用于截断边界x = d 上就意味着没有反射波( 垤方向传播的波) 。于 是，方程( 2 2 0 a ) 可以看成是截断边瓢= d 上的理想的吸收边界条件。 对v 作t a y l o r 展开，有： c 1 一三c 岳+ 和+ ， 亿2 2 ， 若仅取第一项近似，则得到m u r 一阶吸收边界条件： ( 0 。一c a 。) fl 。；o = 0 ( 2 2 3 ) 若取前两项近似，则可得m u r 二阶吸收边界条件： ) ) ) 9 a b 1 o o 乞 2 2 ( 2 2，、，l 南京邮电大学硕上学位论文 第二章f d t d 算法理论 ( a ；a ，一c o t 2 + ( a ；+ a ：) ) fb = 0 ( 2 2 4 ) 将m a r 吸收边界条件中的偏导数用差商近似就可得出对应的时域差分方程。以场分 量为例，m u r 一阶、二阶吸收边界条件对应的差分公式如下： e ；n + l ( 。，露+ 圭) = 霹( ，七+ 三) + c a t - a x 彰“( - ，尼+ 1 一彰( 。，露+ 圭) ( 2 2 5 ) 一e ：t ( o ，七+ 吉) = 一e ( 1 ，| | + 专) + 竺筹尝 彰+ 1 ( 1 ，后+ ) + 彰一1 ( 0 ，后+ ) + 盖 霹( 吖肿 ) + 霹( 1 _ 七十 ) +鹄霹(。，歹+1，尼+吉)一2e(。，歹，足+吉)+霹(。，j一1，霓+吉) + e ( 1 ，j + l ，七+ 圭) 一2 e ( 1 ，j ，尼+ ) + e ( 1 ，一1 ，尼+ ) ( 2 2 6 ) + 霹( o ，豇+ 号) 一2 e ( o ，七+ 圭) + e ( o ，j f ，k 一号) + 霹( 1 ，尼+ 号) 一2 彰( 1 ，尼+ 吉) + 霹( 1 ，k 一 ) 二阶以上的边界条件很少用。 2 2 算例 2 2 1 均匀平面波对介质平面垂直入射的f d t d 模拟 设一均匀t e m 平面波沿z 方向传播，当电磁波在传播过程中遇到不同媒质时，必然 有一部分电磁能量被反射回来，而另一部分将能穿越分界面透射出去，本算例将通过 f d t d 法模拟这一过程。 电磁场量在传播过程中与x 、y 无关。即昙：o ，晏：o 。这是一个一维情况。于是 m a x w e l l 方程变为： 一警= 占拿o t + 啦 ( 2 2 7 ) 一 2 占_ + 础j ，一、 d z ( 2 2 7 ) 誓2 叫警一吗 在第四十个空间网格处加一个正弦平面波s i n ( 2 幸p i 幸f f ) ，f 取为1 0 0 m h ：。空间步 长a - - 0 3 。时间步按照a t = a ( 2 * e ) t 玟值。c 为在自由空间中波传播速度。第五十空间步 南京邮电大学硕士学位论文 第二章f d t d 算法理论 之后为一无耗介质，即仃( m ) = o ，脾= l ，q = 8 5 。1 0 0 时间步的模拟结果如图2 2 ： 刚 忖 ；vvv 】 v 01 0锄加 5 0 7 n 即即t 图2 2 均匀平面波垂直射入介质平面的幅值分布( 1 0 0 时间步， 横坐标0 - - , 1 0 0 代表空间网格，纵坐标代表场幅值) 由图2 2 可以看出波向左已经传播到达边界，由于边界处设置了一阶m u r 吸收边界条 件，几乎没有反射，介质边界处，在介质左边由于介质的反射，波源( z = 4 0 ) 右边的波 幅增大，波源左边的波幅减小，而介质右边的透射波波幅、波速都有减小，这非常符合实 际情况。 图2 3 均匀平面波垂直射入介质平面的幅值分布( 1 0 0 0 0 时间步， 横坐标0 1 0 0 代表空间网格，纵坐标代表场幅值) 继续让波前进，我们给出1 0 0 0 0 个时间步的结果，见图2 3 。 可以看出随时间的增加，介质边界左边的自由空间中的波幅度已经大于波源的最大幅 值1 ，接近于1 5 ，右边介质中的波幅相对于介质边界总体上变小，但是波幅有一个周期 1 4 南京邮电大学硕：七学位论文第二章f d t d 算法理论 性的变化，最大值已经接近于l 。这个幅值是截断边界和介质边界反射叠加的结果。我们 给出相同计算环境下不设介质边界计算1 0 0 0 0 时间步的结果，如图2 4 ，发现没有介质边 界时具有良好的吸收特性，这说明图2 3 的幅值变大是由于介质边界的设