（信号与信息处理专业论文）时域有限差分法及其在电磁波测井正演问题中应用.pdf

摘要 摘要 电磁波测井正演问题指：在已知地层模型电参数和激励源情况下，用 电磁场数值计算方法模拟复杂介质中电磁波的传播，进而定量计算出电磁 波测井仪器的测井响应，本质是复杂介质中电磁场的计算问题。电磁波测 井正演是测并资料反演的重要基础，正演结果的分析也可为仪器参数优化 设计提供重要依据。时域有限差分算法( f d t d ) 是一种行之有效的电磁场数 值计算方法。 由于电磁波测井地层模型是横向的地层，同时要模拟的电磁波测井仪 器是圆柱形，所以本文主要研究圆柱坐标系下，复杂介质中f d t d 法及其在 电磁波测井响应正演模拟中的应用。本文分析并解决了应用f d t d 法的相关 技术性难题包括：数值计算的稳定性条件；f d t d 算法的激励源的分类和设 置的改进；边界吸收条件的处理及设置等。运用f d t d 算法最大技术性难题 是如何设置强有力的吸收边界条件，以保证计算的等效性和精度，吸收边 界条件的吸收效果和吸收精度有待于进一步提高。本文研究了f d t d 算法的 非分裂场完全匹配层( u p m l ) 司题，详细给出完全匹配层的图示及内部电磁 参数的设置，推导出非分裂场完全匹配层的f d t d 表达式，并数值实验结果 证明了u p m l 层的吸收效果明显优于m u r 吸收层。推导出适合求解电磁波测 井响应的u p m l f d t d 表达式，与已有的f d t d 法在电磁波测井正演模拟中 的应用相比较，理论上分析了改进完全匹配层的时域有限差分法的优越性。 本文在柱坐标系下f d t d 法的研究及求解方案的建立，为电磁波测井响 应模拟做好充分前期准备。 关键字电磁波测井；时域有限差分法；吸收边界条件；完全匹配层；正演 模拟 燕山大学工学硕士学位论文 a b s t r a c t f o r w a r ds i m u l a t i o no fe l e c t r o m a g n e t i cw e l l l o g g i n gt o o l si nc o m p l e x b o r e h o l ee n v i r o n m e n t si su s i n go fe c c o n 廿i ec o m p u t a t i o n a lm e t h o dt os i m u l a t e t h ew a v eo fe l e c t r o m a g n e t i ct r a v e l i n ga tk n o w np a r a m e t e r so f m o d e ls t r a t u m a t t h es a n l et i m e ，w ec a n c o m p u t a t et h er e s p o n d i n g o fe l e c t r o m a g n e t i c w e l l l o g g i n gt o o l s i nf a c ti ti se l e c t r o m a g n e t i cf i e l de s t i m a t i o np r o b l e m s t h e s i m u l a t i o no fe l e c t r o m a g n e t i cw e l l - l o g g i n gt o o l si sai m p o r t a n tf o u n d a t i o no f t h er e c o v e r , t h et o o l so p t i m i z e dd e s i g n i n ga r ea l s ob a s e do nt h es i m u l a t i o n r e s o l u t i o n s f d t d ( f 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 ) h a sb e e np r o v e nt o b eah i g h l ye f f i c i e n tt e c h n i q u e s i n c et h es t r a t u mm o d e lo fe l e c t r o m a g n e t i cw e l l l o g g i n gi st r a n s v e r s e ，a t t h es a n l et i m e ，t h ee l e c t r o m a g n e t i cl w dt o o li sa l s oc y l i n d r i c a l ，s ow es t u d yf d t d i nc y l i n d r i c a lc o o r d i n a t e o n eo ft h e g r e a t e s tc h a l l e n g e so fa p p l y i n gt h i s t e c h n i q u ei st oo p e nr a d i a t i o np r o b l e m s ，t h i sh a sb e e nt h ed e v e l o p m e n to f a c c u r a t ea n dc o m p u t a t i o n a l l ye f f i c i e n ta b s o r b i n gb o u n d a r yc o n d i t i o n s e x i s t i n g a n a l y t i c a la b s o r b i n gb o u n d a r yc o n d i t i o n ( a b c ) ，s u c ha sm u ra n dg e d n e y s p e r f e c t l ym a t c h e dl a y e r 口m l ) i nt h i sp a p e r , n u m e r i c a le x p e r i m e n t s a r e a d d r e s s e df o rc o m p a r i n gt w ok i n d so fa b s o r b i n gb o u n d a r yc o n d i t i o n ss u c ha s m u ra n du p m l ，t h o s ec o m p a r i s o n ss h o wt h a tt h eu p m l t e c h n i q u ew o r k sb e t t e r t h a nm u ra b c a na n i s o t r o p i c - m e d i u m ( u u s p l i t ) p e r f e c t l ym a t c h e dl a y e r ( u p m l ) f o r m u l a t i o ni sp r e s e n t e df o rt h et r u n c a t i o no ff d t dl a t t i c e sw i t h g e n e r a ll o o s e l ya n dd i s p e r s i v em e d i a t h en u m e r i c a li m p l e m e n t a t i o no ft h e p m li s s t r a i g h t f o r w a r d w i t h o u t a n ys p l i t t i n g o ft h e f i e l d s ，l o w e r i n g c o m p u t a t i o n a lc o s t s b ya n u m e r i c a l e x a m p l e ，i ti s s h o w e dt h a tu p m l p e r f o r m a n c ei ss u p e rt h a nt h eg e n e r a lp e r f e c t l ym a t c h e dl a y e r t h e n , w ei n f e r o u tt h ee x p r e s s i o no fu p m l - f d t di nc y l i n d r i c a lc o o r d i n a t e s i nt h e o r y , w e n a b s t r a e t a n a l y z et h es u p e r i o r i t yo f u p m l f d t d t h o s es t u d y i n gw o r k so ff d t di n c y l i n d r i c a l c o o r d i n a t ea r e f u l l y p r e p a r i n gf o rf o r w a r ds i m u l a t i o no f e l w da tt h ee a r l i e rp e r i o d k e y w o r d se l e c t r o m a g n e t i cw e l l l o g g i n g ；f i n i t e d i f f e r e n c et i m ed o m a i n ( f d t d ) ；a b s o r b i n gb 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 ； f o r w a r ds i m u l a t i o n i i i 燕山大学硕士学位论文原创性声明 本人郑重声明：此处所提交的硕士学位论文时域有限差分法及其在 电磁波测井正演问题中应用，是本人在导师指导下，在燕山大学攻读硕士 学位期问独立进行研究工作所取得的成果。据本人所知，论文中除已注明 部分外不包含他人已发表或撰写过的研究成果。对本文的研究工作做出重 要贡献的个人和集体，均已在文中以明确方式注明。本声明的法律结果将 完全由本人承担。 作者签字- , v ，j ，2 -日期：阳萨牛月修日 燕山大学硕士学位论文使用授权书 时域有限差分法及其在电磁波测井正演问题中应用系本人在燕山 大学攻读硕士学位期间在导师指导下完成的硕士学位论文。本论文的研究 成果归燕山大学所有，本人如需发表将署名燕山大学为第一完成单位及相 关人员。本人完全了解燕山大学关于保存、使用学位论文的规定，同意学 校保留并向有关部门送交论文的复印件和电子版本，允许论文被查阅和借 阅。本人授权燕山大学，可以采用影印、缩印或其他复制手段保存论文， 可以公布论文的全部或部分内容。 保密口，在年解密后适用本授权书。 本学位论文属于 不保密d ( 请在以上相应方框内打“”) 作者签名：马日纱 导师签名。彳勃芦 导师签名：删7 琴l 。乃， 日期：知叮年午月冶日 日期：如呷年争月冶日 第1 章绪论 1 1 课题的研究背景 第1 章绪论 电磁波测井( e l e c t r o m a g n e t i cw a v el o g g i n g ，e m w l ) i l l 作为地球物理测 并重要手段之一，能够同时测量井周围地层的电导率和介电常数，而这些 常数对于评价地层的含油性是至关重要的，普通电阻率测井和感应测井都 是利用地层孔隙流体的导电性质来区分含油气和含水地层的。电磁波测井 的源频率从几百千赫到几十兆赫，一般采用单发双收探头测两接收线圈的 感应电动势之间的相位差西和幅度比s 。这种相对量，受井眼影响小且空 间分辨率高，广泛地应用于复杂油气地层的评价和随钻测井( l o g g i n gw h i l e d r i l l i n g ，l w d ) 的地质导向中。电磁波测井可分为几兆赫兹的传播电阻率 测井1 2 1 和几十兆赫兹的高频电磁波测井 3 1 ，传播电阻率测井只提供地层电 阻率信息，而几十兆赫兹的高频电磁波测井能同时提供地层介电常数和电 阻率两者的信息，在判别水油层上独具优势，同时也有极好的识别薄层的 能力，高频电磁波测井使用频率要比传播电阻率测井高得多，其纵向分辨 能力比传播电阻率测井的好，但其致命的弱点是它的探测深度不如传播电 阻率测井的深。 几兆赫的传播电阻率测井和几十兆赫的高频电磁波测井的测井仪器结 构和测量原理是相同的，如图1 1 所示，图中的非对称仪器由一个发射线圈 瘌两个接收线圈r l 和恐组成，发射线圈通常加以源频率为f = 2 n z o 的电 流，电磁波测井仪器居于井轴上，测量两个接收线圈感应电动势的相位差 舻和幅度比s 。按习惯，将r 1 和r 2 的中点作为测量记录点，该点到发射源 z 韵距离上称为仪器的源距，r 1 和r 2 之间的距离工称为仪器间距，描述电磁 波测井仪器性能的参数是( 厂，l ，工) 。图1 1 中的对称化仪器有两个对称放 置的发射线圈兀和死，测井时乃和死交替发射，把它们交替发射时的两次测 井响应进行平均后的结果作为测量值，因此对称化测量与非对称化的测量 燕山大学工学硕士学位论文 相比，它有利于减少井孔不规则对测井响应造成的影响，但测井仪器变长 了，测井响应的正演计算也相应增加了一倍，增加了资料处理和反演的时 间。 tr lr 2 t i r lr 2 t 2 ll 图1 - 1 对称和非对称电磁波测井仪器示意图 f i g 1 - 1t h es y m m e t r i ca n du n s y m m e m cc o n f i g u r a t i o no f t h ee l e c t r o m a g n e t i cl w d t o o l s 电磁波测井研究的主要问题是正反演的问题，反演问题是电磁波测井 资料处理本质的问题。电磁波测井反演问题一般可表达成如下式。 y = ，( x ) + e( 1 - 1 ) 记个测井响应为y 维实数域响应空间中的矢量j ，记个地层参数为麟实 数域空间中矢量工。电磁波测井资料处理本质的问题成为从雕参数空间到 艘匡响应空间的逆映射，电磁波测井的正演计算可以看成从艇匡参数空间到 】，维响应空间的映射。正演的意义在于，电磁波测井正演的计算结果可对 仪器参数设置进行优化改进，更好地提高仪器的性能。同时，数值模拟的 正演结果决定着电磁波测井资料反演的水平，是电磁波测井资料反演的基 础和核心依据。本文主要求解电磁波测井的正演问题中的正演方法。 1 2电磁波测井的正演问题和正演方法 电磁波测井正演问题是指在己知地层参数和激励源的情况下，用电磁 场数值计算方法模拟电磁波在介质中的传播过程，电磁场数值计算的方法 2 第l 章绪论 很多，包括有限元法( f i n i t ee l e m e n tm e t h o d s ，f e m ) 、频域有限差分发( f i n i t e d i f f e r e n c ef r e q u e n c yd o m a i n ，f d f d ) 、传输线矩阵法( t r a n s m i s s i o nl i n e m a t r i x ，t l m ) 、数值模拟匹配层法( n u m e r i c a lm o d em a t c h i n g ，n m m ) 及傅 立叶一汉克共轭梯度变化法( c o n j u g a t e g r a d i e n tf a s tf o u r i e r - h a n k e l t r a n s f o r mt y p em e t h o d s ，c g f f h t m ) 、时域有限差分法( f i n i t e d i f f e r e n c e t i m e - d o m a i n m e t h o d ，f d t d ) 等，都己广泛应用到电磁波测井正演模拟中。 1 2 1 正演问题 电磁波测井环境下地层的非均质性 3 5 1 表现如图【2 】1 2 所示。实际地层 电磁参数分布情况是非常复杂的，由于地层沉积过程中形成的纵向不均匀 性和钻井泥浆侵入造成的径向不均匀性，往往呈现地层横向多层不均匀的 分布，甚至可以出现宏观各项异性的分布。电磁波测井反演问题中，电磁 波测井资料的解释除依赖于反演方法外，还要借助于电磁波测井的正演结 果。电磁波测井正演是测井资料反演的基础，决定着电磁波测井资料反演 的水平。斜井中电磁波测井仪器测井模型如图1 - 2 ，劝发射线圈，发射线 圈加入激励源的情况下，用电磁场数值计算方法模拟电磁波的传播，计算 出两个接受线圈电场或磁场值，进而求出代表电磁波测井响应的幅度比和 相位差，这就是电磁波测井的正演过程。 图1 _ 2 倾斜井径中电磁波测井仪器测井模型 f i g 1 - 2c o n f i g u r a t i o no f t h ee l e c t r o m a g n e t i cw e l l l o g g i n gt o o li n s i d ei n c l i n eb o r e h o l e 电磁波测井正演问题就是指：在已知地层模型电参数和激励源的情况 燕山大学工学硕士学位论文 下，用电磁场数值计算方法模拟复杂介质中电磁波的传播，进而定量计算 出电磁波测井仪器的测井响应，本质上是复杂介质中电磁场的数值计算问 题。本文研究的是时域有限差分法应用于电磁波测井响应正演模拟中，用 f d t d 法模拟电磁波在复杂介质中传播，进而计算出电磁波测井仪器两个 接受线圈电动势的幅度比和相位差。 1 2 2 正演方法 目前为止，许多学者都对电磁波感应测井和电磁波测井的数值模拟方 法进行了大量研究，取得了很有意义的成果。水平井中阵列感应成像仪的 响应和倾斜的层情况下深感应和中感应的测井响应【6 ，1 7 】研究已很多。m a r k 等计算了各种复杂地层中感应测井的响应。g r a c i e t 和b a r b 盯a 等人利用有限 差分( f i m t ed i f f e r e n c e ，f d ) 的方法，c h e w 和d a v i d 分别利用了水平特征模 式和垂直特征模式匹配方法计算了单个水平界面感应测井响应。这些方法 都属于频域计算方法，一次只能计算单一频域的数据。频域有限差分法和 有限元法都是直接离散m a ) ( w e l l 方程，有限元法作为一种数值计算方法已 经在工程技术的许多领域内得到了广泛的应用，它可以较方便地解决工程 中许多过去感到棘手甚至无法解决的问题，已成为求解固体力学、热传导 问题、流体流动问题和电磁场、声场等问题的一种有效的计算方法。有限 元法可在圆柱坐标系中模拟测井仪器，然而往往需要求解大型的稀疏矩阵， 给求解带来困难。数值模拟匹配层算法、傅立叶汉克共轭梯度变化澍8 9 】 具有高效快速等优点被广泛应用电磁波测井正演模拟中，但是它们不能应 用于任意测井仪器中。时域有限差分法是一种发展了的有限差分方法，它 将传统的空间差分加入对时间的差分，从而对m a ) 【w e l l 方程进行时间和空 间直接离散，随时间推进交替计算电场和磁场值，f d t d 具有较直观的物 理意义解释。本文讨论柱坐标下时域差分法在电磁场计算中的可行性，并 比较了几种不同吸收层模型的吸收效果。时域有限差分法作为一种数值算 法，对电磁波测井正演模拟有重大的意义。时域有限差分法叶氏( y e e ) 首先 提出的，近些年在电磁波测井正演模拟中的应用【m b 】受到国内外研究人员 的广泛重视。它在电磁波计算中的普遍应用性源于它自身的一些优点如下： 4 第1 苹绪论 不受物理模型几何形状的限制；不受制于散射因素；是一种对瞬变问题的 时间域直接解法；便于理解易于掌握。本文推导出柱坐标系下时域有限差 分表达式，验证了不同匹配层吸收效果，初步建立了电磁波测井正演模拟 的f d t d 法。在己知地层参数和激励源的情况下，用f d t d 法模拟电磁波在 复杂解之中传播过程，本文模拟电磁波测井仪器的测井过程，通过分析计 算结果定量分析出测井仪器的测井原理，为电磁波测井资料的反演奠定坚 实的基础。 1 3 论文研究的主要内容和结构安排 本论文内容主要是研究适合电磁波测井正演模拟的时域有限差分算 法，研究时域有限差分法的技术性问题如：激励源、稳定性条件、吸收层 的问题等。用时域有限差分法模拟电磁波在复杂介质中的传播过程，进而 计算出电磁波测井响应。本文研究重点是推导出柱坐标系下u p m l f d t d 算法的表达式，解决了时域有限差分法吸收层的问题，实验比较了完全匹 配层吸收效果的优越性，为电磁波测井正演模拟做好充分的前期准备。 全文共分五章，具体的章节内容安排如下。 第1 章绪论介绍论文研究的背景，概述了当前电磁波测井的意义及测井 仪器的发展，阐述了模拟电磁波测井正演模拟数值计算方法及意义，最后 介绍了本文的主要研究内容和论文的组织。 第2 章讨论了关于时域有限差分法的一些基本知识，包括：时域有限差 分法的发展、现状、特点。推导了时域有限差分法基本方程，对时域有限 差分法相关技术做了相应要求( 稳定性、数值色散特性) 。 第3 章总结时域有限差分法中激励源类型和设置问题，推导出磁偶极子 激励源的电磁场的表达式。 第4 章研究了吸收边界条件及其处理方法。为推广时域有限差分法应用 范围着重研究吸收边界条件问题，包括吸收边界条件类型如m u r 吸收边界 条件、p m l 吸收边界条件、u p m l 吸收边界条件。重点讨论了圆柱坐标系 下的p m l 层的时域有限差分格式。 5 燕山大学工学硕士学位论文 第5 章详尽推导了适合电磁波测井的圆柱坐标系的u p m l f d t d 法的 表达式，通过编程验证了u p m l 层吸收效果的优越性。在模拟电磁波测井 响应中，初步建立了电磁波测井正演问题的f d t d 求解方案。 最后，对整个论文工作做了总结，并探讨了下一步需要完善和进一步 研究的问题。 6 第2 章时域有限差分法 第2 章时域有限差分法 2 1时域有限差分法的简介 时域有限差分法【l 扣1 ”是求解电磁问题的一种数值技术，它是在1 9 6 6 年 g q k s y e e 第一次提出4 j 的。f d t d 法直接将有限差分式代替麦克斯韦时域 场旋度方程中的微分式，得到关于电磁场分量的有限差分式，用具有相同 电参量的空间网格去模拟被研究体，选取合适的场初始值和计算空间的边 界条件，可以得到包括时间变量的麦克斯韦方程的四维数值解，通过傅里 叶变换可求三维空间频域解。 2 1 1 时域有限差分( f d t d ) 算法的简要发展和现状 近些年来，在电磁计算及电磁工程的各个领域，迅速发展起来了一种 直接求解m a x w e l l 方程组的时域有限差分算法，该算法是对未知的近区电 磁场进行时间和空间上的采样，用有限差分式代替m a x w e l l 旋度方程的微 分式，选取合适的场初值和计算边界吸收条件，得到m a x w e l l 方程的时域 解。k s y e e 首次提出了一种合理的电磁场空间网格划分方法和著名的y e e 氏网格，建立了时域有限差分法。该方法将时间和空间离散化，并将电场 和磁场在时间和空间上交替地分布在各个网格单元上，每个电场分量的周 围都由磁场分量环绕着，每个磁场分量周围都有电场分量环绕着，这正好 符合 m a x w e l l 方程组的旋度方程( 法拉第电磁感应定律和安培环流定律的 微分形式) 所体现出的物理意义，再用中心差分形式近似代替m a x w e l l 方程 组中的微分形式，将m a x w e l l 方程组中的两个独立旋度方程化成差分方程。 通过对网格节点的电磁参数及一定的初始条件的设置，并在网格空间中反 复迭代运算差分方程，可以模拟电磁波传播及电磁波与目标的相互作用过 程，获得m a x w e l l 方程组时域解，并通过傅里叶变换可以获得频域解。在 k s y e e 建立时域有限差分法的基本框架之后，经过许多学者不断完善，时 7 燕山大学工学硕士学位论文 域差分法逐渐走向成熟。尽管时域有限差分法是在研究电磁散射问题时被 提出来的，但在八十年代以后又被广泛地应用到电磁波散射问题研究之外 的电磁兼容、微波电路设计、天线辐射、生物电磁学及电磁波测井等电磁 场工程的各个领域。 随着f d t d 被人们广泛地接受并应用，f d t d 方法本身在不断地完善和 发展，主要表现在以下几个方面。 ( 1 ) 总场边界条件的引入由于引入总场边界条件，使得f d t d 计算区 域划分为总场区域和散射场区。这样的好处是任意模拟波的入射变得简单 易行，无论是平面波或瞬态波，使得只吸收单向波的吸收边界条件设置得 以实现，方便远场的外推计算。 ( 2 ) 远近场的变换f d t d 的模拟只能限于有限空间，为能获得计算域 以外的散射场，必须借助等效原理应用计算区域以内的近场散射数据实现 计算区域与以外的远场的外推。 ( 3 ) n 格剖分网格剖分技术的应用提高了f d t d 方法的有效性和准确 性。这些新技术的应用大大拓宽了f d t d 方法的应用范围。 ( 4 ) 吸收边界条件吸收效果的提高为了用有限的计算区域模拟无限 自由空间的电磁问题，必须在计算边界上设置吸收边界条件，将计算区域 截断。吸收边界从开始简单的差值边界，到后来广泛采用的m u r 吸收边界 条件，近几年发展起来的完全匹配层吸收边界，吸收效果越来越好。1 9 9 4 年，b e r e n g e r 提出完全匹配层吸收边界条件，完全匹配层由于与以往的边 界相比吸收效果有了很大的提高，在f d t d 的计算中被广泛采用。是近几 年吸收边界领域研究的热点问题。 2 1 2 时域有限差分法的特点 随着高速电子数字计算机的发展，计算电磁学也迅速地发展起来，各 种复杂电磁问题的解决提供了有力的手段。在诸多电磁场的数值计算方法 中，时域有限差分法是近年来发展最为迅猛、最受广泛关注和应用范围最 广的一种方法他具有很多突出的特点。 时域有限差分法自提出，在以后经过众多学者的不懈的努力，使之不 第2 章时域有限差分法 断的完善，现已基本成熟。但是在许多方面仍在继续发展，解决问题的能 力和应用范围仍在不断地提高和扩大。 时域有限差分法有以下几个方面的特点。 ( 1 ) 直接时域计算时域有限差分法直接把含时间变量地m a x w e l l 旋度 方程在y e e 氏网格空间中转化成为差分方程。在这种差分格式中每个网格 点上的电场( 或磁场) 分量仅与它相邻的磁场( 或电场) 分量及上一时间步该 点的磁场( 或电场) 有关，在每一时间步计算网格空间各点的电场和磁场分 量，随时间步的推进，就能模拟电磁波的传播及其与物体相互作用的过程。 时域有限差分法把各类问题都作为初值问题来处理，使得电磁波的时域特 性被直接反映出来。应用这一特点使它能直接给出非常丰富的电磁场问题 的时域信息，给复杂的物理过程描绘出清晰的物理图像。如果需要频域信 息，则只需对时域信息进行f o u r i e r 变换。为获得宽频的信息，只需在宽频 谱的脉冲激励下进行一次计算。 ( 2 ) 广泛地适应性在时域有限差分法的差分格式中，被模拟空间电磁 性质的参量是按空间网格给出的。因此，只能设定相应空间点以适当的参 数，就可以模拟各种复杂的电磁结构。媒质非均匀性、各项异性、色散特 性和非线性等均能很容易的进行精确模拟。 ( 3 ) 节约储存空间和计算空间在时域有限差分法中每个网格电场和 磁场的六个分量及其上一步的值是必须储存的，此外还有描述个网格电磁 性质的参数及吸收边界条件和连接条件的有关参量，它们一般与空间网格 总数成正比。由于每个网格的电磁场都按同样的差分格式来计算，所以， 所需要的主要计算时间而言，也是与成正比。相比之下，用矩量法计算 时所需要的存储空间与( 3 加2 成正比，所需的c p u 时间也与( 3 帅胛成正比， 而对f d t d 法而言，两者均只与成正比，所以，节省了大量的存储空间 和计算空间。 “) 适合并行计算以直角坐标系中立方体网格空间为例，若每个坐标 方向的网格数为力个，则计算网格的空间总数n = n 3 。若用厅珂6 个处理器， 则每个处理器只需记忆和处理一行中一个场分量的有关信息，聆斫亍同时 处理，对于一个确定的时间步，全部运行时间就正比于完成一行处理器所 9 燕山大学工学硕士学位论文 需时间，这个时间又正比于一行中一个场分量的个数胛。可见，并行计算可 使的时域有限差分法所需的储存空间和计算时间都成倍的减少，提高了 f d t d 法解决实际复杂问题的能力。 ( 5 ) 简单、直观、易于掌握时域有限差分法直接从m a x w e l l 旋度方程 出发，不需要任何导出方程，简便易懂，在时域中模拟出电磁波的传播及 其与物体作用的物理过程。这种方法容易得到推广，并在很广泛的领域发 挥作用。 2 2f d t d 法的基本方程 在电磁场数值分析的计算方法中，有限差分法是应用最早的一种方法， 它具有简单、直观的特点，一直得到广泛的应用。运用有限差分的概念求 解问题，其思路是把连续变量分解为离散变量，从而使得偏微分方程求解 转化为差分方程求解。 2 2 1有限差分格式的概念 有限差分法简称差分法基本原理是用差分方程替代微分方程，从而将 复杂的微分方程转化为代数方程，可很方便地应用计算机来进行求解。 设函数如) ，其独立变量有个很小的增量h 定义为a x = h ，相应的该 函数的增量为v = ， + 矗) 一f ( x ) ，它称为函数船) 的一阶差分，因为是两 个变量限量的相差，故称为有限分法。而一阶差分舷) 除以增量h 的值称 为一阶差商。 一阶差商公式如下。 笪：f ( x + h ) - f ( x ) ( 2 1 ) a xh 根据一阶导数的定义，可得公式如下。 八加霉= 牌警( 2 - 2 a x ) r w x 可以看出，当公式( 2 一1 ) 中的h 越小时差分和微分的值将越接近，两者的差 1 0 第2 章时域有限差分法 值将越小。公式( 2 1 ) 称为前向差分，此外还有以下两种差分形式。 后向差分公式为： 一af f ( x ) - f ( x - h ) ( 2 3 ) a xh 中心差分公式为： 笠：f ( x + h ) - f ( x - h ) ( 2 4 ) 缸2 厅 三种差分方法中，采用中心差分的截断误差最小，这可以用泰勒公式 展开来说明，如下所示。 f ( x + h m + 而差+ 刍矗2 睾+ - ( 2 5 ) f ( x - h ) - ，一嗉+ 刍矗2 警+ ( 2 - 6 ) 由公式( 2 5 ) 和( 2 6 ) 可得公式( 2 - 7 ) ，表达式为： f ( x + h ) 一f ( x - h ) _ 2 唾+ 争窘+ 函7 ) 比较公式( 2 7 ) 和公式( 2 4 ) 可以看出用中心差分替代微分两者的截断 误差为o ( h 3 1 ，当h 较小时，用差分代替微分将获得足够的精度。 2 2 2 直角坐标系下时域有限差分格式 在时域有限差分算法中，首先把计算空间划分成网格，来离散空间上 连续分布的物理量，并且只计算网格结点上的物理量。时域有限差分法用 于计算电磁场问题，它基于m a x w e l l 旋度方程，计算电场和磁场的六个分 量，需要多个空间网格结构相互关联。 y e e 在提出f d t d 算法时，给出著名的y e e 氏网格结构图，类似于图2 1 所示的网格体系。y e e 氏网格图很直观的表示出在直角坐标系中网格空间 划分的特点，坐标平面上每个电场分量被磁场环绕，同时每个磁场分量被 电场环绕，电场和磁场分量在空间上交叉放置。这种电磁场的空间结构符 合电磁感应和电波传输的规律，每一个网格单元上的电磁场的结构满足法 燕山大学工学硕士学位论文 拉第感应定律和安培环流定律。这样时域有限差分法就很容易模拟电磁波 的空间传播、辐射及散射体相互作用等物理过程。 图2 - 1y e e 网格单元矢量场分布图 f i g 2 - 1p o s i t i o no f t h ef i e l dv e c t o rc o m p o n e n t sa b o u ta y e eu n i t ec e l l 为了描述计算空间内的各部分介质的电磁性质，需要计算出每个电场 和磁场分量所在空间相应介质的电磁参数。由于m a x w e l l 方程组是支配宏 观电磁现象的一组基本方程，同时时域有限差分算法也是基于m a x w e l l 方 程，它是由微分形式的m a x w e l l 方程出发进行差分离散所得，故推导时域 有限差分法的方程表达式应从麦克斯韦旋度方程出发。 麦克斯韦旋度方程如下所示。 v x 日：竺+ j ( 2 8 ) 研 v x e ：一票一l ( 2 - 9 ) 第2 章时域有限差分法 式中，e 为电场强度( 伏特米) ；d 为电通量密度( 库仑米2 ) ；日为磁场强 度( 安培米) ；口为磁通密度( 韦伯，米2 ) ；崩电流密度( 安培米2 ) ；厶为磁流 密度( 伏特米2 ) 。 一 电磁参数本构关系为： d = e e ，b = a h ，= o e ，厶= 口 ( 2 1 0 ) 式中，s 为介电常数；p 为磁导率；d 为电导率；仃。为磁导率，盯和仃，分 别是介质的电损耗和磁损耗。所以，只要给定每一网格相应媒质的电磁参 数，就使得电磁波对复杂结构的模拟变得方便。 f d t d 算法y e e 网格中电场e 和磁场昭分量空间结点和时间步取值如 下表2 1 所示。 表2 1y e e 网格中e 、h 各分量节点的位置 t a b l e2 - 1t h en o t el o c a t i o no f t h eea n dh c o m p o n e n t sa b o u tay e eu n i t ec e l l 空间分量取样 电磁场分量时间轴诹样 x 坐标y 坐标z 坐标 e x i + 1 2 j k e 节点 e ， j + l 2 k e zj k + z 2 m j + i 2k + 1 2 h 节点 屿 i + 1 2j k + 1 2 n + 1 2 h z i + 1 2j + 1 1 2k 仕且用坐体系甲，化m a x w e l l 旋厦力程展升q j 得死组电能物明力程表 达式如下所示。 鲁一孕0 2 = s 鲁+ 盯乓 ( 2 ) j 一= + f e i ，1 ll 跏讲 警一掣= e 季+ 鸣( 2 - 1 2 ) o x o t磁 。 孚一警=s军4-dxo t 仃e “= = ，r ，1 l w 燕山大学工学坝士掌位论文 等一季叫警乜以- 0 2 “= 一“o f r l z 1 4 - 挪讲 警一等叫警一吒影( 2 - 1 5 ，a z瓠 j a t “y 。 竽一等=一_“孥一o-,g1-10xo t ( 2 - 1 s ) = 一_ 二= 一“- 一 i z l n l 卯 对于公式( 2 - 8 ) 和公式( 2 9 ) 的f d t d 差分离散，中心差分的基本方法是 将空间沿三个坐标上分成很多个网格单元，用缸，缈，止分别表示单元格沿 三个轴向的长度，用缸表示时间增量。在直角坐标系中，网格单元的顶点 坐标( x ，弘z ) ，记为( f ，- ，k ) = ( i a x ，j a y ，k a z ) 。 令f ( x ，y ，：，t ) 代表电场或磁场的某一分量，在时间和空间中的离散取 以下公式为：f ( x ，弘z ，) = 八i a x ，y a y ，k a z ，n a t ) = f ”( f ，j ，k ) 。 用中心差分法近似代替函数f ( x ，弘z ，t ) 在时间和空间一阶偏导数，即 用中心蒡分代替一阶偏导数，公式如下所示。 o f ( x , y , z , t ) 。：坠：! ! 墨! 二：垡二：丝2 缸x o f ( x , y , z , t ) 。：盟：盟二：盟：盟 砂4y o f ( x , y , z , t ) 。：业：生生2 二：业：! 二曼 出业 堑! 兰! ! ! 趔。垫：盟二：垫：丛2 西f ( 2 1 7 ) ( 2 1 8 ) f 2 - 1 9 ) 佗2 0 ) 将公式( 2 1 1 ) 按上述中心差分公式的差分规则分别展开，可得电场值的 表达式如下。 e “( f + ，_ i ) = c a ( m ) e 2 ( ，+ ，七) 十c b ( m ) 墨兰：生：墨! 二墨兰兰垡：二圭：墨! 一 缈 嘭峙( f + ，七+ ) 一巧n q ( f + ，k 一 ) & f 2 2 1 ) 第2 章时域有限差分法 式中： c a ( m ) = ( 1 - 等等) ( 1 + 号等等) ( 2 - 2 2 ) ( 聊) = ( i 赫) ( 1 + 号等筹) ( 2 2 3 ) 将公式( 2 1 4 ) 按上述差分方法，得磁场在工方向的表达式如下。 月：+ ( f ，+ ，_ i + 号) = c 尸( m ) 月：4 ( i ，+ ，量+ ) 一c q ( m ) 式中： 一 t ”z 2 4 ) yi z 叫 e ( f ，_ ，+ ，七+ 1 ) 一髟( f ，+ ，di zj c 尸( 脚) = ( 1 一警) ( 1 + 哥筹) c q ( m ) - _ w _ ，- ) ( 1 + 筹) y ，：方向的e ，日值同理可得。 2 2 3 柱坐标系下时域有限差分格式 佗- 2 5 ) ( 2 2 6 ) 当被研究的对象具有曲线表面时，如天线、光纤、电磁波测井等，若 用直角坐标系下网格单元拟和曲线表面，会形成如图2 2 所示的阶梯形边 界，产生的误差太大。 图2 - 2 梯形元拟合曲线边界 f i g 2 - 2s i m u l a t i n gc u l 、r eb o u n d a r yw i t hs t e pc e l l s 使用这种阶梯形边界会产生两个问题：沿导体介质表面可能激起表面 波传输，引起附加的数值色散；为了拟和曲率半径小或多层媒质表面，就 燕山大学工学硕士学位论文 要减小网格尺寸，相应地时间步长也减小，这既增加了计算存储量，又增 加了计算时间。为了克服阶梯表面所引起的模拟的困难，较好的方法就是 使用曲面坐标系、球坐标系、抛物线坐标系及椭圆坐标系等正交曲线坐标 系，下面将研究圆柱坐标系下f d t d 方程【1 3 】的建立。 圆柱坐标系下的f d t d 网格结构与直角坐标系y e e 网格相同，电场和磁 场分布如图2 - 3 所示。 图2 3 圆柱坐标系网格单元的矢量场分量位置 f i g 2 - 3p o s i t i o no f t h e f i e l dv e c t o rc o m p o n e n t sa b o u tac y l i n d r i c a lu n i tc e l l 圆柱坐标系的坐标变量为( ，驴，z ) ，用于建立柱坐标系f d t d 方程的基本 单元，网格单元的顶点坐标定义如下。 ( ，z ) = ( i a r ，妒，k a z ) - - a t ( ，七)( 2 2 7 ) 式中，a r ( i ) = r ( i + 1 ) 一，( f ) ；( - ，) = 妒( _ ，+ 1 ) 一( ，) ；a z ( k ) = z ( k + 1 ) 一z ( k ) 。 把麦克斯韦旋度方程在柱坐标系中展开，它们等同下面的六组标量方 程组，简便起见只推导出其中三个电场方向的表达式，电场和磁场之间的 关系表达式如下所示。 j a , 。e ：三等一譬( 2 - 2 8 ) ro 口 j o ) e e , ：譬一冬( 2 - 2 9 ) 1 6 第2 苹时域有限差分法 j a , s e ：三掣1 冬( 2 - 3 0 ) ru ll o q , 啦峨：掣一粤( 2 - 3 1 ) ， 巾u h - ：萼一姿( 2 - 3 2 ) - j o ) 1 i h ：一1 了0 ( r e ) 一等( 2 - 3 3 ) l o iu m 利用傅立叶变化的微分特性可将频域场值化到时域场中，傅立叶变化： j j a | 氆q 3 4 ) 取时亥1 t = n + l 2 ，将肠分量转换到时域中。应用中心差分格式法可得 到电场迭代f d t d 方程，电场分量的取值与磁场分量取值在时间和空间上 相互关联，时间和空间取值都是相差半个步长。电场分量的取值公式如下。 五芦o ，+ ，_ i + 曲= f ( _ ，+ ，| j + ) + f s 垡兰! ! ：! ! 墨i ! 二墨蔓! g ：! 生2 一 f 击 目：q ( f ，_ ，+ ，七+ 1 ) 一。日( f ，- ，+ ，砷 业 e r l ( f + ，_ ，七+ ) = 譬( f + ，后+ ) + a t s 。h ，r q 、- 。 x 2 ，女+ 1 ) 一月? q ( f + ，后) & h t + ( i + 1 ，_ ，k + 9 一彬q ( f ，j + ) ， 霹“( f + ，_ ，+ 吉，七) = 霹( f + ，+ ，七) + a t _ s ( f + 1 ) 月：+ 壬( f + 1 ，+ ，| | ) 一日? q ( f ，+ ，_ j ) a r 月：q o + ，j + l ，七) 一日：畸( f + ，歹，_ j ) 击 f 2 - 3 5