（电磁场与微波技术专业论文）fempml算法在电磁散射特性分析中的应用.pdf

西北r 业人学硕十学位论文 摘耍 摘要 有限元方法是计算电磁学的主流方法之一，对复杂结构和非均匀介质问题有 很强的描述能力。在开域电磁散射问题中，有限元需要虚拟边界以截断无限大空 间。本文用完全匹配层作为开域边界条件，对导体目标的敖射问题进行了研究。 论文中首先介绍了结点有限元和矢量有限元的分析步骤，比较了两者的优缺 点，然后引入平面完全匹配层的概念，这部分是论文的理论基础。在应用方面， 研究了平面匹配层矢量有限元法在二维、三维电磁散射分析中的应用，及共形 匹配层结点有限元法在二维散射分析中的应用，分别采用二阶三角形及矢量四 面体作为二维、三维问题的剖分单元。在分析过程中，深入探讨了完全匹配层的 使用条件和参数设置。并通过算例，对比了吸收边界条件与完全匹配层的吸收性 能。文中几个典型散射体的计算结果与其它文献或理论结果的对照表明了这种方 法的有效性和准确性。 关键词：有限元，完全匹配层，电磁散射 塑! ! 三些奎耋堡圭兰堡篁兰 垒! 墼堡垒 a b s t r a c t t h ef i n i t ee l e m e n tm e t h o df f e m ) i so n eo ft h em a j o rn u m e r i c a lm e t h o d si n c o m p u t a t i o n a le l e c t r o m a g n e t i c s d u et o i t s v e r s a t i l i t ya n df l e x i b i l i t y , t h ef e mi s c a p a b l eo fm o d e l i n gc o m p l i c a t e ds t r u c t u r e sa sw e l l a si n h o m o g e n e o u sm a t e r i a l s w h e ns o l v i n ge l e c t r o m a g n e t i c ( e m ) s c a t t e r i n gp r o b l e m s 谢t l lt h ef e m ，t h ei n f i n i t e r e g i o ne x t e r i o rt ot h es c a t t e rm u s tb et r u n c a t e db yaf i c t i t i o u sb o u n d a r y i nt h i sp a p e r , t h ep e r f e c t l ym a t c h e dl a y e r ( p m l ) i sa d o p t e da st h et r u n c a t i n gb o u n d a r yf o rt h e a n a l y s i so fc o n d u c t o re l e c t r o m a g n e t i cs c a t t e r i n gp r o b l e m s i nt h i sp a p e r , i ti n t r o d u c e st h ee s s e n t i a lp r o c e s s e so ft h en o d a lf i n i t ee l e m e n t m e t h o da n dt h ev e c t o rf i n i t ee l e m e n tm e t h o d ，r e s p e c t i v e l y t h et h e o r yo ft h ep l a n e p e r f e c t l ym a t c h e dl a y e r si sa l s ob e e np r e s e n t e d p l a n ep m l v e c t o rf e mi se m p l o y e d f o rt w o a n dt h r e e - d i m e n s i o n a le ms c a t t e r i n g p r o b l e m s ，a n dt h ec o n f o m a a l p m l n o d a lf e mf o rt w o d i m e n s i o n a lp r o b l e m s s e c o n d - o r d e rt r i a n g l en o d a lf i n i t e e l e m e n t sa n dt e t r a h e d r a le d g ef i n i t ee l e m e n t sa r eu s e df o rc o m p u t i n gs c a t t e r i n gf r o m t w o a n dt h r e e d i m e n s i o n a ls t r u c t u r e s ，r e s p e c t i v e l y i nt h i sp a p e r , w ec o m p a r et h e a b s o r b i n ga b i l i t yb e t w e e nt h ea b s o r b i n gb o u n d a r y c o n d i t i o n s ( a b c s ) a n dt h e p e r f e c t l ym a t c h e dl a y e r s ( p m l ) a l s o ，w ed i s c u s sd e e p l y o nt h ea p p l i c a t i o n c o n d i t i o n sa n dt h ep a r a m e t e rs e t t i n go ft h ep m l n u m e r i c a lr e s u l t sa r ei ne x c e l l e n t a g r e e m e n tw i t ht h o s ep r e s e n t e db yo t h e rp a p e r so re x a c ts o l u t i o n s k e y w o r d s ：f i n i t e e l e m e n t m e t h o d ( f e m ) ，p e r f e c t l y m a t c h e dl a y e r ( p m l ) ， e l e c t r o m a g n e t i cs c a t t e r i n g i i 塑! ! ：三些奎兰堡圭i 丝丝塞 堡= 童丝垒 1 1 有限元概述 第一章绪论 有限元方法对复杂结构和非均匀介质问题有很强的描述能力。它建立了一条 连续系统离散逼近的自然途径，这是有限元方法最重要的特性。有限元法的基 本概念是将由偏微分方程表示的连续函数所在的封闭场域化分成有限个小区域， 每一个小区域上用一个选定的近似函数来代替，于是整个场域上的函数被离散 化，由此获得一组近似的代数方程，并联立求解，以获得该场域中函数的近似解。 有限元法和其它求解边值问题近似方法的根本区别在于它的近似性仅限于相对 小的子域中，它将函数定义在简单几何形状的单元域上，不考虑整个定义域的复 杂边界条件，这是有限元法优于其它近似法的一个主要方面。 1 9 4 3 年，c o u r a n t 2 1 首次明确的提出了有限元的思想一最小位能原理与分片 插值的离散形式。1 9 6 0 年，有限元的名称正式出现【3 j 。六十年代中期，数值分析 学家意识到有限元思想的重要性，将偏微分方程理论、泛函分析、逼近论等引入 到有限元的理论体系，建立了有限元方法的数学基础。自此，有限元开始在工程 计算中得到广泛应用，直至成为某些学科数值算法的基础。 1 9 6 9 年，p _ es i l v e s t e r 有关波导模求解的文章是个里程碑，标志着电磁学有 限元方法的出现【4 1 。随后，有限元法被用于求解电机磁场、静电场、波导本征值、 涡流场、散射与辐射场等等，涉及电磁学的各个领域。到九十年代，有限元法己 经同矩量法、时域有限差分法共同构成计算电磁学的三大主流算法。 大量科技文献与商业化软件的出现表明，在电磁学相关领域，有限元理论已 经趋向成熟，然而教射问题的有限元理论是个例外，三十年来这个领域新的方法 层出不穷，但是一直未发展到接近工程应用的水平。显然，这是一个广阔而复杂， 值得深入研究的领域。 1 2 选题背景和意义 在低频区和谐振区，电磁散射的基本分析方法是基于微分方程和积分方程的 数值求解方法【4 9 1 。 西北工业人学硕十学位论文第一章绪论 夺微分方程法 求解散射问题的微分方程数值问题方法主要有两种：一种是有限差分法 ( f d ) ，将连续的三维空间用网格划分开，m a x w e l l 方程化为具有一定精度的差 分方程，在这些网格点的未知电场强度便由代数方程给出。f d 方法分为时域有 限差分( f d t d ) 和频域有限差分( f d f d ) 。另一种方法就是有限元法( f e m ) 。 这两种方法在求解散射问题时，都需要人为地将求解空间限定为有限区域，使未 知量数目保持有限。但是相比而言，有限元方法的计算精度要高于有限差分法， 且有限差分法在模拟复杂几何形状对的误差也要大于有限元法。 夺 积分方程法 在积分方程中未知量是导体表面或涂敷阻抗面的面电流( 面积分方程) 或可 透入散射体内部的体电流( 体积分方程) ，通过应用等效原理体积分可以化为面 积分方程。积分方程的求解域限于散射体表面或内部，离散化后未知数的数目较 微分方程法大为减小。散射场用辐射积分求出，计算精度有保证。最基本的方法 就是矩量法( m o m ) 。但用矩量法得到的代数方程组系数矩阵为稠密的( 即大多 数矩阵元素不为零) 。对于高维矩阵的求逆工作要求有较大的计算机内存且很耗 机时，这是矩量法不能用于大尺寸三维物体的主要障碍。 经过比较，在本文中选择有限元方法来求解电磁散射问题。 有限元方法求解电磁散射问题的困难之处在于，有限元是一种区域性方法， 受计算机内存空间和运行时间的限制；而电磁散射是一个开域问题，这就需要特 殊的边界条件来为有限元确定一个有限大的计算区域。有限元开域边界条件1 5 可以分为局域边界条件( l o c a lb o u n d a r yc o n d i t i o n ) 与全域边界条件( g l o b a l b o u n d a r yc o n d i t i o n ) 。 有限元应用于散射问题之初，研究人员就十分清楚局域与全域边界条件的缺 陷，新的开域边界条件往往给研究工作带来巨大的进展。电磁散射问题的有限元 理论与其它学科相比有许多共同之处，其个性部分都直接或糊接地与开域边界条 件联系在一起，这一点是理解其发展趋势的关键【2 ”。 夺全域边界条件 电磁散射早期的开域边界条件几乎都是全域边界条件，最早引入的是边界积 分方程( b o u n d a r yi n t e r g r a le q u a t i o n ) 。单矩法( u n i m o m e n tm e t h o d ) 是稍后出现 的全域边界条件，在单矩法的边界上，散射场展开为本征函数和的形式，待求变 量是各个奉征函数的展开系数。由于这种方法的边界形状仅限于圆形、球形，故 未得到推广。 有限元矩量法，有限元边界元等方法，在本质上都j 以归入有限元边界积 西北工业大学硕士学位论文 第一章绪论 分法，到八十年代，该方法逐渐得到关注，是公认的精确数值方法。 如前所述，全域边界条件准确地描述了散射体外部空间的电磁特性，是精确 的边界条件，然而由于考虑了边界上所有结点间的相互作用，全域边界条件离散 化之后会形成一个稠密矩阵，使有限元与全域边界条件的耦合矩阵失去稀疏性， 不利于存储和求解。限制了有限元求解电大尺寸散射问题的能力。 夺局域边界条件 吸收边界条件a b c ( a b s o r b i n gb o u n d a r yc o n d i t i o n ) 是典型的局域边界条件， 七十年代末出现，在有限元1 6 、时域有限差分法中得到广泛应用。二维的吸收边 界条件已经有了完整的理论体系1 7 】1 8 1 ，并在实践中广为采用。三维矢量吸收边界 条件目前仍在发展中【9 1 ，有许多问题尚未解决，是研究的前沿方向之一。由于提 高吸收边界条件的阶数或者采用“共形”技术能够减小物体与边界之间的距离， 提高计算效率，所以高阶和共形吸收边界也是目前研究的热点1 1 0 】。现有的高阶和 共形吸收边界条件公式繁琐，计算迭代次数多，未能达到节约计算时间的目的， 稳定性与吸收效果亦不理想，总体上说是不成功的】。 m e i ( m e a s u r e de q u a t i o no f l n v a r i a n c e ) 【50 l 由m e i 等在1 9 9 2 年提出，最初用 于截断有限差分网格，1 9 9 3 年开始被用于有限元开域问题1 5 ”。在m e i 方法中， m e i 系数将边界上某一点的场与其相邻结点处的场联系起来，不会改变有限元矩 阵的稀疏性，且具有能贴近散射体放置的优点。迄今为止，m e i 方法的主要应用 限于二维，公开文献中也没有见到有关稳定性的结论，m e i 的许多优点也需要在 理论上进一步证实。 完全匹配层p m l ( p e r f e c t l ym a t c h e dl a y e r ) 是近年出现的最出色的局域边 界条件。在理想状态下，p m l 对入射波的吸收与频率及入射角无关。不会改变 有限元系数矩阵的稀疏性，而且更好的吸收效果只需要增加p m l 的厚度即可。 大量的研究证明 3 8 埘】，p m l 能够放置在距散射体很近的距离上，且吸收效果优 于a b c 。 完全匹配层的概念由b e r e n g e r 于1 9 9 4 年酋先以分离场变量的形式在f d t d 中提出l ” l 引，这种形式的p m l 对麦克斯韦方程有较大的改动，不能简单的应用 到有限元中。1 9 9 5 年，s a c k s 等提出了直角坐标系中张量形式的p m l 1 4 】。这种 p m l 相对于有限元来说只是一种各向异性媒质，从而得到了广泛的应用。张量 形式的柱面、球面【1 5 _ 埔 p m l 也出现了。最近的研究已将平面、柱面、球面p m l 组合使用【l ”，从而增强了p m l 的共形性，减小了p m l 与散射体问的缓冲区， 并使边缘效应问题有了解决的可能，使计算效率和精度得到了进步的提高。 夺总结 西北上业大学硕士学位论文 第一章绪论 全域边界条件是精确的边界条件，可是破坏了有限元系数矩阵的稀疏性；局 域边界条件只考虑相邻结点间的相互作用，因此保持了系统矩阵的稀疏性，不足 之处是局域边界条件一般为精确辐射边界条件的某种近似，精度和应用范围有 限。另外，局域边界条件必须距散射体一定距离，保持一定的形状，不免扩大了 网格区域；全域边界条件则可以紧贴散射体设置，形状也比较灵活。 局域边界条件的发展方向是缩小网格区域，降低反射误差；全域边界条件的 研究则侧重于有限元系统矩阵的处理和求解。但两者最根本的目标是相同的，即 提高有限元求解电大尺寸散射问题的能力。 在这种情况下，张量形式的完全匹配层在有限元方法中的优越性就很明显 了：第一，可以很容易地应用到有限元的源代码中，并不改变有限元系数矩阵的 稀疏性；第二，可以距散射体尽可能地近，从而可以大大减小计算量，提高计算 效率。特别是第二个优点，可以提高有限元求解电大尺寸散射问题的能力，具有 现实工程意义。 综上所述，为了进一步提高有限元方法求解电大尺寸物体电磁散射问题的能 力，本硕士论文着重研究张量完全匹配层在有限元方法中的实现。 1 3 论文结构安排 本论文的主要研究内容是基于有限元，完全匹配层算法的目标电磁散射特性 分析。根据有限元算法在解决开域散射问题时的方式及完全匹配层的形式，论文 主要内容大致分为四个部分：首先以二维结点有限元算法为例介绍了有限元的基 本步骤；其次介绍了三维结点有限元及矢量有限元算法的基本步骤，在对比中分 析了结点有限元方法在解决三维问题中存在的困难，及矢量有限元的优点；再次 引入了平面完全匹配层的概念，结合三维矢量有限元算法分析目标的电磁散射特 性；最后将匹配层概念扩展到共形匹配层，结合二维结点有限元分析目标的电磁 散射特性。 全文分为五章，主要安排如下： 第一章课题研究的背景、选题意义及国内外的发展状况。 第二章以二维结点有限元算法为例，介绍了有限元分析的基本步骤和思路。 第三章分别介绍了三维结点有限元和矢量有限元算法，在具体的分析过程 中，通过对比，总结了结点有限元方法在解决三维问题中存在的困 难，说明了矢量有限元方法的优点。 第四章引入了张量形式的完全匹配层的概念。通过公式推导，论证了平面 匹配层的吸收机理，说明了匹配层结构参数的构造原则，并通过数 西北= 业大学硕士学位论文 第一章绪论 值例子分析了不同参数对匹配层吸收性能的影响。在此基础上，用 三维矢量有限元，平面匹配层方法来分析目标的电磁散射特性，并与 其它算法得到的结果进行比较。 第五章在第四章引入平面匹配层吸收机理的基础上，讨论了柱面、球面匹 配层的结构参数的构造，分析了不同参数对匹配层吸收性能的影响， 并引申至共形匹配层。在此基础上，用二维结点有限元共形匹配层 方法来分析目标的电磁散射特性，证明了共形匹配层的有效性。 西北工业大学硕士学位论文 第二章二维结点有限元方法介绍 2 1 引言 第二章二维结点有限元方法介绍 目标电磁散射特性的计算是电磁学中的一个复杂问题，普遍采用的策略是先 从二维问题入手，然后再将成熟的方法推广到三维。同时，在实际工作中，若目 标的纵向尺寸远大于横向尺寸，经常用二维模型逼近原来的三维问题。这样做， 一方面可以简化分析、计算过程；更重要的是，可以大幅度的减少对计算机存储 空间的需求和减小计算机的计算量。因此，对二维目标散射特性的分析方法进行 研究，既是研究策略的需要，又有实际应用背景。本文正是基于这样的思路，首 先从二维目标电磁散射特性的研究入手，讨论二维有限元方法。 结点单元在处理三维电磁问题时遇到了困难：但在精度、产生的未知数数目 等方面优于矢量基单元。而另一方面，在应用结点有限元分析二维问题时并没有 遇到伪解问题，因此，在本章中，采用结点单元来处理有限元问题。 2 2 边值问题 根据分析结果j 2 l ，2 2 1 ，当散射体中包含材料或导体的尖角时，在尖角处横向场 具有奇异性。所以，应用有限元方法分析包含材料或导体尖角目标的电磁散射问 题时，需要处理场的奇异性。通常的做法是在了解场的奇异性的基础上选择相应 的特殊插值基函数。但只有几种几何结构的奇异场表达式是已知的；再者，这种 方法难以编制可用于不同问题的通用程序。 然而，在各向同性介质中，任何二维场都可以分解为只有纵向电场的t m 波 和只有纵向磁场的t e 波，而在材料或尖角的纵向场不存在奇异性【2 1 2 2 。因此， 选择以纵向场作为变量的标量亥姆霍兹方程作为场的控制微分方程，可避丌场的 奇异性给有限元带来的困难，从而能够编制处理不同问题的通用程序；这样做还 有一个好处；能避开结点有限元法处理不连续介质分界面的不方便。然而考虑到 有限元方法的普通性，本章以任意的二阶微分方程为例，来描述有限元问题i ”j 。 考虑二阶微分方程所定义的边值问题， 两，| t i j k 大学硕士学位论文 第二章二维结点有限元方法介绍 一丢( q 芸卜昌卜詈 + 舯= ， ( w ) e q c z 2 t ， 式中，巾是未知函数，a x , a 。和p 是与区域物理性质有关的已知参数，f 是源 或激励函数。常用的二维拉普拉斯方程、泊松方程和亥姆霍兹方程是( 2 ，2 1 ) 式的 特殊形式。 这里，所考虑的边界条件为： = p 在r l 上( 2 2 2 ) 以及 【q 罢x 6 + c l y 票多卜徊：g 在r 2 上( 2 2 3 ) 【q i百。y p + 徊2 9 在r ：上( 2 - 2 3 ) 式中，r ( = r 1 + r ：) 表示包围面q 的轮廓或边界，磊是外法向单位矢量，、p 和 q 是与边界物理性质有关的已知参数。特别地，p 和q 能被看作边界源或边界激 励。 假如，q 和a ，表征的区域性质具有不连续性或突变，并且在不连续界面上没 有任何类型的面源，则。满足连续性条件 = 巾一 在1 1 d 上( 2 2 4 ) ( 警三+ 哆等；) 拓一。良a 。- x + 等歹 ：矾上： 式中，f d 表示不连续的界面，上标“+ ”( 或“一) 表示相应的量位于r d “+ ( 或 “- ”侧) ，n 表示r 。的单位法向。( n i n2 1 ) 图2 1 具有不连续界面r d 的区域 哆，p 西北工业大学硕十学位论文第二章二维结点有限元方法介绍 2 3 变分公式 等价于上述边值问题的变分问题可表示为 艇p 叫一乱上 眨s ， 【o =在r 1 上 相应的泛函为 即，= 衅( 黔蚓2 州卜亿。埘 + 北班。中户叫t ，中m j f c 中，= 三f c l ( 罢) ( 警) + q ( 詈) ( 等) + 印砷 q 亿，m + f ，( 归一g ) 删i 一他伪d q 占f ( 。) = 豫卜昙( q 罢) i 导( q 詈) + 肿一十d q + 啦f ( q 罢；) + ( q 割- = ；卜r 仁，川 + f ：( 冲一q ) 鼬d r 州卟蜂一誊豺加：、厂卜亿，固 龇q 罢；) + 一茅多h 归刊陋= 。一 西北丁业大学硕士学位论文第二章二维结点有限元方法介绍 一昙( a 。罢 一号( q 詈 + 肿一，= 。 c z 3 s ， ( q 割+ ( d ，割永阳2 。 b s 刀 显然，它们分别是( 2 2 1 ) 式和( 2 2 3 ) 式。这告诉我们，在这个问题中，( 2 2 2 ) 式是必要的边界条件，它必须用显式方法加上；而( 2 2 3 ) 式是自然边界条件，能 自动且隐含的满足。 2 4 有限元分析 利用有限元方法求解电磁场问题，首先将求解区域分割为有限个单元，然后 分片插值。提高网格模型对原始目标几何模型的逼近程度和场变量模型描述真实 的场变量分布的能力，可以有效地提高计算结果精度。因此单元类型的选择十分 重要。 有限元单元可分为线性单元和高阶单元。高阶单元得到的解比用线性单元所 得解精确的多；一般情况下，单元的阶数越高，解的精度越高。并且，误差近似 正比于。”，这里n 表示未知量的个数，n 表示单元的阶数1 2 4 1 。但是，高阶单 元公式复杂、带宽随阶数增加，线性单元公式相对简单、相应的有限元系数矩阵 带宽窄。综合考虑，在本章二维问题中，选择二阶三角形单元作为基本单元类型。 2 4 1 区域离散和单元插值 1 ) 区域离散 有限元分析的第一步是将面区域q 划分成若干三角形单元，这就是区域的离 散。离散的一个基本要求是单元之间既没有重叠又没有削隔，另外，单元应通过 它们的顶点相连，换句话说，一个单元的顶点只能位于和它相邻单元的顶点上， 而不能位于其它单元的边上。较好的离散应注意以下两点：第一，为减小误差， 应避免窄形单元的生成，尽可能采取等边三角形；第二，对所期望的精度，应保 持单元数最少，比较好的做法是：在解急剧变化的区域应用较小的单元，而在解 变化较平缓的区域应用教大的单元【2 3 j 。 为了识别每个单元，我们用一组整数给单元编码；同样，为了i i 别单元顶点 处的结点，可以用另一组整数给单元编码。因为每个单元与六个结点有关，所1 i ， 一个结点除了有它在整个系统中的位置外，它还有在相应单元中的位置，这种位 堕! ! 三些奎茎堡圭耋堡垒塞 丝二童二丝鱼生童堡垂壅堡坌丝 置也可用整数编码，这种编码称为局部编码。为了将全局结点编码、局部结点编 码和单元编码联系起来，在此引入一个6 m 的整型数组n ( i ，e ) ，其中 i = l ，2 ，3 ，4 ，5 ，6 ；e 一1 , 2 ，3 ，m 。m 表示单元总数。在连续数组n ( i ，e ) 中，i 是结点的局 部编码，e 是单元编码，n ( i ，e ) 的值是结点全局编码。显然编码是不唯一的。另外， 为了便于处理( 2 2 3 ) 式，还引入了类似数组n s ( i ，s ) ，将r ：上的线段和有关结点联 系起来。同时，为了强加狄利克雷边界条件，还需要引入能存储位于r 上结点 全局编码的向量n d ( i ) 。 2 ) 单元差值 在q 域中，第e 个单元内场可以表示成二次函数： 。( x ，y ) = a 。+ 6 。x + c 。y + d 。x 2 + p 。x y + f 。y 2 ( 2 4 1 ) 矿，酽，c 。，d 。，e 。，f 8 为待定的常系数，e 是单元编码。在六个结点上强加( 2 4 1 ) 式，即可确定上式中的待定系数。将它们代回到( 2 4 1 ) 式中，得到： 6 旷x ，y ) = 孵x ，y ) ； ( 2 4 2 ) 卢l 假设结点的编码如图2 2 所示，则上式中的插值函数为 巧( z ，y ) 2 ( 2 巧- 1 ) 弓 ，= 1 ，2 ，一 ( 2 鲫 孵x ，y ) = 4 茸骂，联( x ，y ) = 4 骂骂，联( x ，y ) = 4 骂置 在上式中，e 由下式给出 弓( 础) = 寺( q e + 牛+ c ；y ) j = l ，2 ，3( 2 1 4 4 4 ) 其中 a l ；一一蠼苟；群= 圬一蝣；彳= 巧一 口；= x ；群一蜡x ； 鹾e = 儿e 一_ y ；c ；= 一 ( 2 4 5 ) 口；e = # 虻一并；蟮= y ? 一蛭； = 恐e 一群 肌吉毛 显然，n ；是满足下式的二次函数， f 1 x t ( x ，) = 岛2 o ( 6 f 一噬吖) ( 2 4 6 ) 面积 z2 j f j ( 2 4 7 ) 塑i ! ：；些奎耋堡圭兰堡篁苎 篓三塞三丝丝皇童堡垂查鎏! ! 型 3 图2 2 二阶三角形单元 2 4 2 里兹方法的计算公式 1 单元公式推导 为了简单起见，我们首先考虑奇次诺曼边晃条件，即( 2 2 3 ) 式在，= q = o 时 的特例。在此条件下，( 2 3 2 ) 式给出的线积分为零。因此，泛函可写为 吖 f ( 中) = f 。( 巾。) ( 2 4 8 ) = 1 式中，m 表示单元总数，f 。是子泛函 叩) = 圳q ( 爿+ 文等) 2 州。叮卜亿。 一心。，蚴q 其中，q 。表示第e 个单元的区域。将( 2 4 2 ) 式带入，并将f 。对吖求导，得到： 弦 p _ ，等1 r ：吒 细；i l + j 1 1 也 其元素如下： ( 2 4 1 0 ) ， 阵篓噬撇竺科 r。l j e1 提 p e 1 j堡彬卅 阵 l 斑 哞 一丌 式 单 耍j ! 三些奎兰堡主兰堡丝塞星三量三堡垒皇查堡垂窒鎏坌塑 k ：fc 】，卜警警+ a y 等v j , i 等+ m 叫哟 b 4 m ， 巧= 他，卜蔷吉万亏邓m 叫哟 ( 2 4 1 1 ) f ，- ，= 1 ，2 ，3 ，6 列向量 b 。 的元素为 彤= m y g ；d x d y f _ 1 ，2 州3 一，6 ( 2 4 1 2 ) 如果在每个单元内的系数a 。、a ，、和源f 均为常数，并分别用a e 、a ；、。和f 。表 示，那么可以用解析法求解( 2 4 1 1 ) 式和( 2 4 1 2 ) 式，在这个过程中需要用到的公 式为 腿。( 研) 协) “( 孵) ”蚴。揣z 岔 得到 弦 = + p 群= 0i = 1 ，2 ，3 式中， 和 伊 是如下对称矩阵 ( 2 4 1 5 ) b = 等f e i = 4 , 5 6 ( 2 4 1 6 ) 巧= 1 4 8 j o 五- - 1 l 咚, q e 叶, + 茚c ) f ，= 1 ，2 ，3 或 或= 髭= 1 3 a 。 1 3 。 l 3 。 耳蛭+ ) 暂巧+ q e q e c 3 e ) e d 2 e 口3 e + q e c 2 e c ；) 筏= 0 以= 爿盖= 爿品 或= 筏，= 砭，畿= 或 以= 去 ( 娣2 + 耳鹾+ ) + ( c ? 2 + 彳+ 2 ) 髭= 去h e e e + 2 b ；b ；+ 娣蟹+ 2 ) + a ；( c ；c ；+ 2 c ；c 3 h ；c 2 ；2 ) 鬣= 去 ( w 蟮+ 2 霹西+ 圩酲+ 耳2 ) + y i _ ，c l 。c 3 。+ 2 c 2 。c 3 。十c i + c p ) 鸽= 去 ( 蟛2 + 鬈+ 蟛) + 日；( 2 + c 2 。c ，e + 2 ) 1 2 里! ! 三些查兰堡主耋堡篁塞 篁三塞三丝丝皇童堡垂壅鎏坌型 畿= 去 ( 蟮w + 2 醛娣+ 甥+ 骘2 ) + ( c 并+ 2 c ；c ? + c ；蝎e 2 ) 心= 去 n ：( 骘2 + 苗钟+ + “胎2 + “埘2 ) 憎 ：坐 la 1 8 0 6一l一1 一l6- 1 一l一16 oo-4 40o o_4o o一40 00_4 _ 4 00 3 21 61 6 1 63 21 6 1 61 63 2 2 方程组的组合 将所有m 个单元组合起来，然后再对f 应用驻点条件，得到方程组： 黔封井辩俐一邛， 亿。加， 可以写成紧凑形式 【k 】 中) = b ( 2 4 1 8 ) 式中， i x = 缸i 。 ， 6 = 势) ( 2 4 1 9 ) 这个组合的过程是将巧加到k 幢咖( 。) 上，将鲜加到瓦( 上。 3 边界条件的迭加 首先考虑第三类边界条件。上面的方程组是在中位于r ，上满足奇次边界条 件的假设下导出的。现在我们考虑( 2 2 3 ) 式中非零y 和非零q 的一般形式。 在这种情况下，( 2 4 8 ) 式的泛函增加一附加项： 吒( 巾) = f ：( 考0 2 - q o ) d r l ( 2 t 4 2 0 ) 假设1 1 ：由m 。个单元边或线段组成，则( 2 4 2 0 ) 式可写成 e ( 中) = 巧( m 5 ) ( 2 4 2 1 ) 类似于一维情况，每一段内的未知函数中可近似为 = j o j f 2 4 2 2 ) 西北工业犬学硕士学位论文第二章二维结点有限元方法介绍 式中 n i = 1 一，n ；= ( 2 4 2 3 ) 其中f 是沿线段从结点1 量到结a2 的归一化距离。具体的蜕，在结a1 ， f2 0 ；在结点2 ，f = 1 ；在两个结点之间，它按线性规律变化。9 4 ( 2 4 2 2 ) a 带 入胃，并将它对巾；求导，得到 鏊= 缸l y 嘲谜一l q 唧( 2 4 2 4 ) 式中，1 5 表示线段的长度。上式也可写成矩阵形式 澍_ k 妒 4 筋， k s 和 中的元素为 k 5 s i 。st l + 0 ) ( 2 4 2 6 ) 鲜= 9 3 ( 2 4 2 7 ) t j 将只考虑进方程组，( 2 4 1 7 ) 式应修正为 料兰o o , 固l j + 射井辩俐一 + 篓( i 5 万) 一 否5 ) = 。 ( 2 4 2 8 ) 为了扩展向量和矩阵，显然需要一个能将线段和有关结点的全局编码联系起来的 数组。这个数组用n s ( i ，s ) ( i = 1 , 2 ；s = 1 ，2 ，3 ，m ；) 表示，它与前面讨论的联系数组n ( i ，e ) 类似，即，用n s ( i ，s ) 存储第s 线段上第i 个结点的全局编码。有了这个数组后， 将每个蟛加到j t n 咖刖一上，即可将 k5 组合进f k 】；同样，将矿加到k 亿，) 上， 即可将 b 5 组合进 6 ) 。 其次讨论狄利克雷边界条件的强加，置大数法是常用的一种简洁有效的方 法。假设结点3 位于r ，上，并有巾，= 见，将这个条件加到方程组中。选择一个 非常大的数，例如1 0 ”，令 k 3 3 = 1 0 ”，b 3 = p ，x 1 0 7 0( 24 2 9 ) 西北工业大学硕士学位论文第二章二维结点有限元方法介绍 则与啦确关的方程变为 k 3 i m i + k 3 2 中2 + 1 0 7 0 m 3 + k 3 4 中4 + 茂5 m 5 + k 3 6 中6 = p 3x 1 0 7 0( 2 4 3 0 ) 假设所有的矩阵元素和未知量均远小于1 0 ”，则( 2 4 3 0 ) 式实际上等效于巾，= p ，。 这样方程组仍保留了对称性。对于r 上具有m 个结点的般问题，当 i = 1 ，2 ，3 ，n i 时，只要 墨d ( ) 州( ，) = 1 0 ”，( 。) = p ( ，) 1 0 7 0 ( 2 4 3 1 ) 就可以加上边界条件。 2 5 本章小结 本章以二维问题为例，简要介绍了标量有限元方法的原理和建立过程，采用 里兹法描述了一般二阶微分方程的有限元分析过程，在此基础上，可以推导出任 意特殊边值问题的有限元解。所选择的单元为二阶三角形单元，是实际二维问题 中常常采用的典型单元。 西北工业人学硕士学位论文第三章三维有限元分析 3 1 引言 第三章三维有限元分析 传统的有限元方法在使用定义在结点上的标量基函数描述三维矢量场时，遇 到了困难。这些困难在一定程度上限制了有限元方法在三维电磁场问题中的应 用。 新近出现的矢量有限元比较好地解决了上述困难。与结点有限元不同，矢量 有限元是使用定义在棱边上的矢量基函数描述电场或磁场。其实，早在1 9 5 7 年 w h i t n e y l 2 s l 就描述过这种类型的单元，但它们在电磁学中的应用和重要性直到近 些年才被认识到。 虽然结点有限元方法存在着不足，为了便于理解三维问题中的有限元方法的 应用，并进一步理解结点有限元方法在解决实际问题中存在的困难，本章首先进 行结点有限元分析。在此基础上，引入了矢量有限元的概念，从而可以通过对比 更明确的了解矢量有限元法在解决三维问题中的优势。 3 2 三维结点有限元分析 在本节中，为了便于说明有限元方法的原理，首先从普通边值问题入手，进 行分析【2 ”。 3 2 1 边值问题 三维边值问题可用如下二阶微分方程定义： 一昙卜票 - 参一；硒- j _ 昙( 吐警 + 肿= 厂一瓦l q i 厂万【q一瓦l 吐i j + 肿叫 相应的边界条件为 中= p 存s 上 ( 3 2 2 ) 西北r 业大学硕士学位论文第三章三维有限元分析 f q 娑6 + a y 票多+ a z x + a y t 细；1 _ = ；+ 椰：g 在马上 ( 3 2 3 ) 【q i面，z 2p + 舯2 9 酗2 上 ( 3 2 3 ( 等，+ o 咖0 + y q 警；卜 ：f 等；+ 丐等多+ 警0 在昌上 。2 5 2 【q i _ 。+ q 百厂y + 哆i 。p 在昌上 3 2 2 变分公式 式中 上面定义的边值问题等价于下列变分问题 f 8 f ( q ) ) = 0 【中= p在s 上 ( 3 2 6 ) ，c 中，= 三m 卜( 罢 2 + q ( 爹) 2 + q ( 罢 2 + 肿 矿 。，2 ， + 压睁2 巾声堋例y 无论a 。、a ，、a ：、和，是复数还是实数，上式均成立。如果这些参数都是实数， 上式可化为 ，c 中，= 吾m q 剖2 + q 矧2 + 呸嘲2 制吖 y + 1 i ( r 1 i 2q m - q c 1 ) p s ( 3 2 8 ) 一言f ( 加+ 伊+ 砂 西北工业大学硕十学位论文第三章三维有限元分析 类似于一维和二维情形，如果d ，、a v d a ：存在不连续性，则( 3 2 6 ) 式应该和连续 性条件( 3 2 4 ) 式一起使用。 3 2 3 有限元分析 三维问题的有限元公式的建立非常类似于二维情形。这里使用线性插值的基 本四面体单元来描述。 1 区域离散 2 图3 1 线性四面体单元 首先将体积v 离散，即将它划分成许多小体积单元，在此应用如图3 1 所示 的四面体单元离散计算域，单元结点局部编码如图3 1 。结果，体积表面s 被划 分成许多小三角形单元。然后，用一组整数给所有体单元编号，用另一组整数给 四面体顶点处的所有结点编号。与二维情况相似，引入整型数组m f j ，e 1 将单元编 码和结点编码联系起来，i = 1 ，2 ，3 ，4 ，e = 1 ，2 ，3 ，m ，m 表示体单元总数。为 了便于强加狄利克雷边界条件，引入整型数组n d ( i ) 来存储s ，面上结点的全局结 点编码。为了处理第三类边界条件，引入3 x m ，的整型数组，用来表示s ，面上的 面三角形单元及与它们相关结点的关系。这个数组用瑚f i ，s 1 表示，这里i = 1 ，2 ，3 ， 而s = 1 ，2 ，3 ，肘。，其中，m 。足s ，上面三角形单元的总数。其它需要的数据包 括：每一结点的坐标每一体单元的q 、口，a z 、掰u f 的值，s 上每一结点的中的 给定值，以及上每一面三角形单元的r 年i j q 值。 2 单元插值 将区域离散后，小阴面体单i 内的未知函数能够近似为 ( x ，y ，z ) = d 。+ 矿x + c 。y + d 。z ( 3 2 ，9 ) 在单元的四个结点上强加( 3 _ 2 9 ) 式，则可以确定上式中的四个系数。将第j 个结 点上的币值记为o ；，则可得 巾：= d 。+ 6 。+ c 。并+ d 。z ： 中；= 。+ 扫。十c 。m e + z ； 拈专匿 k 其中 拈击 c 。= 嘉 d e ：上 6 v 。 o ；= 矗8 + 6 。b e + c 。蝣+ d 。 m ：= 口。+ 6 2 + c 。成+ d 。荔 o ；巾；m ； 蚝 兵螺 巧e z ： 矿t = 三 6 戎 中； = 嘉( 叩e e 一“+ 嚆咿巾；) 击( 槲+ + 鹕+ 蝌) = 击( 州+ 琏西；+ 鹕e e + “) = 嘉( 和i + 磷蝼+ 咿瞒) 。单元体积 将行列式展开后，可以确定系数。；、垮c ，e i ，e 。 将、扩、c 。和d 。的表达式代回( 3 2 ，9 ) 式，得到 ，1 9 喵虻兹， 噬蝣， 哦虻， 哦，瞒，噬， 薪吖彳 。乃喵 虻噶 片吖 蛾 ，菇 蠼，茸k彳 ! 塑! ! 三些丝堡圭兰堡篁塞 量三塞三丝童堡重坌堑 ( t y ，z ) = n ，e ( x ，y ，z ) 中： ( 3 2 1 0 ) 式中插值函数j ( x ，y ，z ) 为 蟛( 圳，z ) 2 矿1 ( 嘭+ 彤x + 丐_ y + 争) ( 3 2 1 1 ) 类似于二维情形，可以证明，插傻函数具有如下性质： n ；( x j , y j , z j ) = 岛= o竺： ( 3 z 1 2 ) 另外，当观察点位于四面体单元的第j 个结点的对面上时，；( 工，y ，z ) 为零。因 此，四面体单元各个表面上的位函数由该面上三个顶点的值决定，从而，单元之 间的连续性得到了保证。 3 里兹方法的计算公式 在离散和插值后，可用里兹方法将问题表示成未知量为结点场值的方程组形 式。 ( 1 ) 单元方程的推导 首先考虑( 2 4 3 ) 式在，= q = 0 的情形下的特例。这种情形下泛函只包括体积 分，可写成 f ( o ) = ( 审。) ( 3 2 1 3 ) 式中m 表示体单元总数，f 。由下式给出 叫叫= 丢i t 卜( 詈) 2 + 嘭( 詈) 2 + 呸( 罢) 2 + ( 叫2 卜。2 。， 一m 。，咖 其中，p e 表示第e 个单元的体积。将( 3 2 1 0 ) 式带入到( 3 2 1 4 ) 式中，然后取对 的偏导数，得到 面o f e2 骞审；肌卜警警乜可o n 可o n ；坦警警+ n f n ) v b 。嘲 一l 妣n ：d v 上式也可写成矩阵形式 ( 3 2 1 6 ) 一中k 啦 j | o 丝渺 ，【 西北工业人学硕十学位论文第三章三维有限元分析 = 甄剑：d v ( 3 2 1 8 ) 显然， 是对称矩阵。如果系数q 、日，、珂：、和勃耖在每一单元内均为常数， 并分别用a ，e 、a e ，、。_ 乖咿。表示，则可以用解析法求( 3 2 1 7 ) 和( 3 2 1 8 ) 式的积 分。在这个过程中，用到一个基本公式： 皿( 研) 。( 婀) 协) ”( ；) ”咖2 高等s v 8 ( 3 2 1 9 ) 结果为 巧= 互专( 鲜+ c ；+ 彬) + 而v e 。( 1 + 吒) ( 3 2 2 。) 筇2 j ，。 ( 3 2 2 1 ( 2 ) 组合成方程组 有了( 3 2 ，1 6 ) 式中给出的单元方程，能够组合所有的m 个单元，然后再对f 强加驻点条件，z 1 4 方程组 = 芒t 丽- r l j = 善啊阡_ o ( 3 z 忽) 其紧凑形式为 【k 】 o ) = 6 ) ( 3 2 2 3 ) 其中，【足】和 6 ) 分别由 i 。 和 否0 组合而成，即将每一个k 加到k n 吐哪一上， 将茸加到。1 上。 ( 3 ) 考虑第三类边界条件 在，o 或可0 的情况下，( 3 2 1 3 )