（通信与信息系统专业论文）基于三棱柱网格的三维矢量有限元方法研究.pdf

摘要基于三棱柱网格的三维矢量有限元方法研究 摘要 有限元方法( f e m - t h ef i n i t ee l e m e n tm e t h o d ) ，是在电磁学研究中应用广泛的一种 计算机辅助分析方法。这种方法用许多子域来代表整个连续区域。在子域中，未知函 数用带有未知系数的简单插值函数来表示，于是，整个系统的解用有限数目的未知系 数近似。然后，用里兹变分或伽辽金方法得到方程组。最后，通过求解方程组得到边 值问题的解。 本文的研究对象是基于三棱柱剖分单元的矢量有限元方法。有限元方法的常用剖 分网格是四面体单元。实际上，三棱柱单元是另一种有用的有限元三维剖分网格，它 对于电场分量在棱柱的延展方向上变化不大的情况，以及建模分层及深腔等结构的目 标特别方便。因为这种剖分网格允许在不规则的二维平面上进行自由剖分，而在第三 个维度上进行规则的剖分。即首先在目标单层或横截面上离散成二维的三角形网格， 然后沿目标的纵向展成三维的棱柱网格。 本文对基于三棱柱剖分单元的矢量有限元方法进行了算法研究。首先，深入研究 了有限元方法的基本原理。在此基本分析框架下，重点研究了当剖分网格为三棱柱单 元情况下的算法实现。最核心工作是单元矢量基函数的给出和单元矩阵的生成。 在算法研究的基础上，本文搭建了基于三棱柱剖分单元矢量有限元方法的分析平 台，包括前处理过程、系统矩阵的生成过程和后处理过程。最后，将该平台应用于波 导腔谐振频率和波导中含有最佳匹配层( p m l ) 问题的数值分析，讨论其表现。 关键词：电磁场，矢量有限元方法，三棱柱剖分单元，波导腔谐振频率，波导中含 有最佳匹配层( p m l ) 问题，数值分析 a b s t r a c t 基于三棱柱网格的三维矢量有限元方法研究 a b s t r a c t n u m e r i c a l t e c h n i q u e s h a v eb e c o m ea m a j o rf o r c e i nt h ef i e l do f a p p l i e d e l e c t r o m a g n e t i c o n eo ft h ei m p o r t a n tt e c h n i q u e su s e di st h ef i n i t ee l e m e n tm e t h o d ( f e m ) i ti sv e r yp o p u l a ri nm e c h a n i c a la n dc i v i le n g i n e e r i n gb e c a u s ei ta l l o w st h e a n a l y s i so fc o m p l e xs t r u c t u r e s i ta l s oi sc o m p u t a t i o n a l l ye f f i c i e n tb e c a u s ei ty i e l d ss p a r s e m a t r i c e s i nt h i sp a p e r ，t h et r i a n g u l a rp r i s me l e m e n t sa r eu s e df o rd i s c r e t i z i n gt h ec o m p u t a t i o n a l d o m a i nw i t he d g e - b a s e dv e c t o rf i n i t ee l e m e n tm e t h o d t od a t e ，t h et e t r a h e d r o ni so f t e n t h ee l e m e n tf o rd i s c r e t i z i n gt h ec o m p u t a t i o n a ld o m a i n i nf a c t ，t h et r i a n g u l a rp r i s me l e m e n t i sa n o t h e rc h o i c ei nn u m e r i c a ls o l u t i o n so fe l e c t r o m a g n e t i cf i e l dp r o b l e m s t h i sk i n do f e l e m e n ti sv e r ys u i t a b l ef o rt h ed i s c r e t i z a t i o no fm u l t i l a y e r e ds t r u c t u r e sa n dd e e pc a v i t i e s ， w h i c ha r eu n s t r u c t u r e di nt w od i m e n s i o n sa n ds t r u c t u r e di nt h et h i r dd i m e n s i o n w e p e r f o r md i s c r e t i z a t i o ni nt h ew a yt h a tt h et r i a n g u l a rp r i s me l e m e n t se x t r u d ea l o n gt h e l a y e r - g r o w t hd i r e c t i o nw h i l ec a p t u r i n gt h ei r r e g u l a rg e o m e t r yi nt h et r a n s v e r s ec r o s s s e c t i o n t h ee d g e b a s e dv e c t o rf i n i t ee l e m e n tm e t h o du s i n gt h et r i a n g u l a rp r i s ma s t h e e l e m e n tf o rd i s c r e t i z a t i o ni sp r o p o s e di nt h i sp a p e r a tt h ev e r yb e g i n n i n g ，t h eb a s i ct h e o r y o ff e mi ss t u d i e d t h eb a s i ca n a l y t i c a lr o u t i n ef o rt h ee d g e - b a s e dv e c t o rf e mu s i n ga n y o ft h ee l e m e n t si ss a m e ，w h i c hi so nt h eb a s i so fc o n v e n t i o n a lf e m h o w e v e r , h o wt o p e r f o r mt h ea n a l y s i so ft h ef e mc o n c r e t e l yi so fg r e a td i f f e r e n c ew h e nd i f f e r e n te l e m e n t s a r eu s e df o rd i s c r e t i z a t i o n t h a ti st os a y , w h a tw ed oi st os t u d yt h ec o n c r e t ew a yt o p e r f o r mt h ea n a l y s i so ff e ms t e pb ys t e p a n dt h em o s ti m p o r t a n tt a s ki st h ed e r i v a t i o no f t h ev e c t o rb a s i sf u n c t i o n sf o rt h et r i a n g u l a rp r i s me l e m e n ta n dt h eg e n e r a t i o no ft h em a t r i x s y s t e m a p l a t f o r mf o rn u m e r i c a ls i m u l a t i o ni sb u i l tb a s e do nt h ea l g o r i t h mw h a tw ep r e s e n t ， i n c l u d i n gt h ep r e p r o c e s s i n gr o u t i n e s ，t h em a t r i xg e n e r a t i o n a n dt h ep o s t - p r o c e s s i n g r o u t i n e s t h e nw ee x a m i n et h ev a l i d i t yo ft h ep r e s e n t e de d g e b a s e df u n c t i o n sb y c o n s i d e r i n gt h ee i g e n v a l u e so ft w od i f f e r e n tc a v i t i e su s i n gt h et r i a n g u l a rp r i s ma s t h e t e s s e l l a t i o ne l e m e n t ，a n dd i s c u s st h ep e r f o r m a n c eo ft h ea l g o r i t h ma n dt h ep l a t f o r m k e yw o r d s ：e l e c t r o m a g n e t i cf i e l d s ，t h ef i n i t ee l e m e n tm e t h o d ( f e m ) ，t h et r i a n g u l a r p r i s me l e m e n t s ，t h ec a v i t yp r o b l e m ，t h ew a v e g u i d ep r o b l e m ，n u m e r i c a la n a l y s i s 声明尸l 刃 本学位论文是我在导师的指导下取得的研究成果，尽我所知，在 本学位论文中，除了加以标注和致谢的部分外，不包含其他人已经发 表或公布过的研究成果，也不包含我为获得任何教育机构的学位或学 历而使用过的材料。与我一同工作的同事对本学位论文做出的贡献均 已在论文中作了明确的说明。 研究生签名：年月日 学位论文使用授权声明 南京理工大学有权保存本学位论文的电子和纸质文档，可以借阅 或上网公布本学位论文的部分或全部内容，可以向有关部门或机构送 交并授权其保存、借阅或上网公布本学位论文的部分或全部内容。对 于保密论文，按保密的有关规定和程序处理。 研究生签名：年月日 硕j ：学位论文 皋十二三棱柱州格的二三维矢量有限元方法研究 1 绪论 1 1 研究背景 使用计算机数值仿真来解决科学及工程问题是在二战后随着电子计算机的发展 而开始出现的。计算机仿真技术，就是首先根据所研究问题的控制方程建立相应的数 学模型，在模型确定之后，选择相应的算法并在计算机上实现。计算机仿真最早的应 用见于1 9 5 2 年，当时被用于天气预报并取得了很大的成功。从此以后，计算机数值 算法得到了迅猛的发展，出现了越来越多的数值仿真方法并涉及到各种不同的领域， 例如流体动力学、物理、化学、天文、电磁学等等。伴随着数值仿真技术的发展产生 了许多新的学科方向。计算电磁学就是以电磁理论为基础，以高性能计算机为辅助工 具，结合高效的数值算法而诞生的- - f q 解决复杂电磁场与微波工程问题的新学科。随 着科学的发展，数值仿真技术在工程应用中占有越来越重要的地位，并且不断有新的 数值算法产生。在当今社会，数值算法已经成为开发新产品、新设备，解决工程问题、 推动科学发展与社会进步的重要手段，它同理论与实验共同构成了当代科学研究的支 柱。 人类的计算能力取决于计算工具的性能与计算方法的性能两个方面。一种高效的 数值算法离不开高性能计算机的支持。如今高性能电子计算机的发展日新月异，就个 人微机而言，2 0 g h z 以上c p u 与1 g 以上内存的微机已随处可见，完全能够满足数 值仿真试验的需要，并且价格二般来说仅需数千元，这为数值仿真技术的应用提供了 很大的便利。然而对于提高计算能力来说，数值算法的计算效率的提高比硬件性能的 提高更为重要。因此自从上世纪5 0 年代以来一直到目前为止，计算方法的研究都非 常活跃，并出现了一系列的数值仿真方法。主要包括偏微分类方法如有限元法、有限 差分法；积分类方法如矩量法、边界元法；以及其它如模匹配方法、射线类方法等等。 各种方法都具有自身的特点和局限性。其中有限元方法由于其强大的模拟任意几何结 构及任意特性的复杂材料的能力，并且具有适合并行处理、程序通用性强等许多优越 的特性，从而获得了人们的喜爱和广泛的应用。另一方面，作为差分类方法，有限元 法导出的为高度稀疏的线性系统。而对于以积分方程为基础的数值方法产生的离散方 程，其系数矩阵通常为稠密矩阵，所有元素都需要大量的数值计算。因此虽然积分类 方法的未知量可能比有限元方法的少，但其数值计算工作量很大。这也是有限元法与 矩量法等方法相比的优势所在。目前有限元方法已经几乎被运用到电磁场中的各个方 i 绪论 硕上学位论文 面，成为数值算法中最重要最成熟的方法之一，并且有许多国内外的学者致力于有限 元算法的研究。近几年来，每年世界上都有数十篇关于有限元研究的科技论文出版， 并且有系统介绍关于有限元技术及在电磁场与微波技术领域应用的书籍，如金建铭的 电磁场有限元方法等。目前许多国家都在有限元方法的研究上投入了很大的人力 物力，并研发出了许多实用成熟的商用软件，如h f s s ，a n s y s ，f e k o 等等，已经 投入广泛的使用，解决了许多的工程问题。 1 2 研究工作的应用意义 我们在运用有限元方法来分析电磁场问题的时候，常用的剖分网格是矩形块单元 和四面体单元。其中，矩形块单元的优点在于它们比较简单，但是主要缺点是只适用 于有限类的几何结构；而四面体单元的优势恰恰就在于对模型几何形状的高度适应 性，它几乎可以处理任意不规则形状的模型。因此，我们常常选择四面体单元来进行 三维网格剖分。 实际上，三棱柱单元是另一种有用的结构剖分网格。因为对于一些具有特定特征 的模型结构来说，这种剖分网格兼顾了运用四面体单元剖分时所带来的形状适应性和 运用矩形块单元剖分时所带来的便捷。具体来说，三棱柱单元允许在不规则的二维平 面上进行自由剖分，而在第三个维度上则进行规则的剖分。也就是说，我们首先在二 维平面上用三角形面元进行剖分，对一个具体的三维模型来说，就是选定该模型的一 个边界表面，用三角形面元进行剖分：然后，完成体剖分的方法，就是以二维平面为 起点以一定的高度成倍延展相同的三角形网格。 那么，上面所说的一些具有特定特征的模型结构又是什么样的呢? 也就是说，有 哪些结构是最适合用三棱柱单元进行离散，并进一步展开分析的呢? 三棱柱单元作为一种有限元三维剖分网格，对于建模分层及深腔结构的目标特别 方便。对于分层结构来说，可以首先在目标单层上离散成二维的三角形网格，然后沿 着层次的延伸方向展成三维的棱柱网格。对于深腔结构来说，则是首先在目标的横截 面上离散成二维的三角形网格，然后沿着目标的纵向展成三维的棱柱网格。 但是，至今，基于三棱柱单元的矢量基函数并没有得到广泛的应用，而由此展开 的理论研究，相比于矩形块单元、四面体单元，以及六面体单元，也显得有些不足。 基于这样的现状，本文对基于三棱柱单元的矢量有限元方法以及应用展开研究，是具 有较强的理论和应用意义的。 2 硕上学位论文基于三棱柱网格的三维矢量有限元方法研究 1 3 本文的主要工作和章节架构 本文研究的矢量有限元方法所基于的剖分网格是三棱柱单元。对于不同的剖分单 元，有限元方法的总体分析流程是相同的，区别之处在于各步骤的实现方法。所以， 基于三棱柱单元的矢量有限元法的研究方法，是在研究矢量有限元方法的基本框架 中，针对三棱柱单元，首先研究每个步骤的具体实现，然后实现各步骤的整合。 下图是矢量有限元法的基本分析流程，以及各步骤中的具体工作内容。 任务 i 工作内容 i2 、对单元 图1 3 1 矢量有限元法的基本分析流程 不同剖分单元的矢量有限元方法都是在上面的基本分析框架中展开的，然而，不 同的剖分单元在算法上又有其各自的特殊性。也就是说，具体到上面各个流程中的每 个步骤，不同剖分单元的实现方法不尽相同。所以，本文的写作是按照从一般到特殊 的思路展开的。下面对本文的章节架构予以说明。 第一章：绪论。从有限元方法的研究背景说起，着重探讨了本文的研究对象 基于三棱柱网格的三维矢量有限元方法的应用意义，概述其研究现状。 l 绪论硕士学位论文 第二章：有限元法的基本原理。以最常用的有限单元四面体单元为例，介绍 了矢量有限元法的基本原理。总的来说，矢量有限元方法是将单元棱边上的切向场分 量作为待求量，用矢量函数构造基函数，直接将场量进行展开，进而求解矢量场。由 于未知量取在棱边上，这种方法也叫棱边元法。 第三章：基于三棱柱网格的三维矢量有限元方法。基于上图所示有限元方法基本 分析框架，研究剖分单元为三棱柱单元情况下矢量有限元方法的实现算法。 在本算法中，针对三棱柱单元，最为核心的研究工作是建立单元s 矩阵和t 矩 阵。包括：建立单元局部坐标系，选择单元矢量基函数，生成s 矩阵和t 矩阵的矩 阵元素。所以，在章节写作中，本文将此核心工作放在首要位置予以阐述。 除了单元s 矩阵和t 矩阵生成的算法研究，当我们编制计算机程序来实现用三 棱柱单元进行区域离散，进而对待求全域展开矢量有限元的分析时，不同的剖分单元 具体的实现方法有别。所以，继系统矩阵生成的算法研究之后，本文将针对基于三棱 柱网格的三维矢量有限元方法中的前处理过程的特殊性予以阐释。 第四章：基于三棱柱网格的三维矢量有限元方法在本征值问题中的应用。 在电磁学中，一类重要的问题是本征值问题，经常出现的波导腔体谐振频率问题 以及在波导中含有p m l 波的传输问题。第四章将基于三棱柱网格的三维矢量有限元 方法应用到上述问题的求解中。数值仿真的目的有二：一是验证算法研究的正确性， 二是验证有限元程序的正确性。 本章从两个方面着手研究。一方面，在算法上，不同的应用问题中矢量有限元方 法的算法有别，具体体现在：单元矩阵的形成，系统矩阵的形成，以及后处理过程。 对此，本文将结合具体的应用问题在第四章中予以阐释。另一方面，在实现上，本章 给出了一些数值算例，验证了本文算法对实际问题的仿真计算能力。 第五章：总结与展望。首先，总结了本文的主要工作。因为本文的主要工作意义 在于，为基于三棱柱网格的三维矢量有限元方法研究搭建最初的工作平台，所以，在 之后的篇幅中给出了可以进步深入拓展的方向。 4 硕i ：学位论文基于三棱柱网格的三维矢量有限元方法研究 2 有限元法的基本原理 2 1 引言 有限元方法( f e m ) 是以变分原理和加权余量法为基础的数值计算方法。其思想最 早出现于上世纪4 0 年代，在5 0 年代被用于飞机设计，然而其开创性的工作是 r w c l o u g h 于1 9 6 0 年完成的，随后这种方法得到快速的发展并广泛地应用于工程中 的结构分析问题中。将有限元法移植到电磁工程领域是二十世纪七十年代的事情。有 限元方法在电磁领域中的应用始于对封闭的电磁系统的分析，这些问题都可以归结为 在给定边界条件下的电磁场边值问题。应用g a l e r k i n 加权余量法导出有限元线性系统 的方法拓宽了有限元法的应用范围，使其适用于具有复杂边界形状与复杂边界条件、 含有复杂媒质的问题。这些年来，有限元法的研究日益深入，与有限元法相关的一些 数值仿真技术也有了很大的进展，使有限元技术能够解决越来越多的电磁问题，如三 维场建模求解、开边界问题、散射问题的求解、白适应网格划分、高磁性材料及具有 磁滞饱和非线性特性介质的处理等，使有限元技术有了很大的发展，并已广泛的应用 于工程问题的仿真，成为解决工程和科学问题的一种通用方法。 在计算电磁学领域，有限元方法的发展经历了从基于结点的有限元法到基于棱边 的矢量有限元法再到高阶矢量有限元法等一系列阶段。其中结点有限元法适合于静电 问题，用于泊松方程中的标量电势的求解，但用于求解高频电磁问题中的矢量电场或 矢量磁场时则会出现一些问题，即不能够保证各单元相邻表面之问场的连续性，并且 不能够正确地表示场的旋度的零空间，从而在有限元仿真过程中有伪解出现，使仿真 结果不可靠。因此，在基于结点有限元法的基础上发展了矢量有限元方法。矢量有限 元方法一般用于求解基于电场或基于磁场的矢量h e l m h o l t z 方程。根据所求解问题控 制方程的不同所采用的基函数也不相同。当用于模拟基于电场的矢量h e l m h o l t z 方程 时，根据电场的物理特性，矢量有限元法采用切向矢量基函数来展开电场，保证了各 单元相邻表面之间电场的切向连续性，而对场的法向连续性不作要求。对于大多数的 电磁场边值问题，一般采用基于电场的矢量h e l m h o l t z 方程进行分析求解。 总的来说，到目前为止在电磁场中应用的有限元方法共有三种类型。即结点有限 元、切向连续矢量有限元、法向连续矢量有限元。矢量有限元克服了传统的结点有限 元的不足，在分析高频电磁问题时取得了很大的成功，因此在电磁场数值计算中与切 向矢量有限元相关的高效数值算法是目前研究的热点方向。 2 有限兀法的幕本原理硕士学位论文 2 2 有限元法的基本原理 2 2 1 电磁场边值问题 用有限元法分析电磁场问题，首先要对电磁场中的边值问题进行研究。边值问题 出现在物理系统的数学模型中，如何对其快速而有效的求解一直是数值模拟方法研究 的主题。典型的边值问题可用区域q 内的控制微分方程和包围区域q 的边界r 上的边 界条件来定义。对于大多数的电磁问题，都能用如下的控制方程来表示： l o = g( 2 2 1 ) 这里三表示微分或积分算子，g 是已知的激励函数，是需要求解的未知场函数， 可以是电场、磁场或电势、磁势等。对于一般的问题，主要有d i r i c h l e t 边界条件， n e w m a n 边界条件、及辐射边界条件等。 计算电磁学所分析的时谐电磁场的特性可以表述为如下的电场波动h e l m h o l t z 方程： v x ( j 1 v 雹( ，) ) 一碍。雹( ，) = 一心z 0 t 7 ( 2 2 2 ) 或磁场波动h e l m h o l t z 方程： v x ( 。v 豆( r ) ) 一瑶以宜( ，) = v c 。歹 ( 2 2 3 ) 对于无源微波器件，上式中歹= 0 。( 2 2 2 ) 与( 2 2 3 ) f l 口为高频电磁问题的控制方程。 一般采用( 2 2 2 ) 式来进行有限元建模。 在时谐有限元分析中用到最多的边界条件为理想电壁边界条件与理想磁壁边界 条件。在理想导电壁表面上，边界条件可以表示如下： 而e ( r ) = 0 o r r f( 2 2 4 ) 在理想磁壁表面上，有： 而v e ( r ) = 0 o n r ( 2 2 5 ) 确定了控制方程与边界条件，就可以使用有限元方法对所分析的问题区域进行离 散，进而求出待求的未知场量及其他物理特性。 2 2 2 伽辽金加权余量法与里兹变分法 因为在电磁学中大多数重要的实际工程问题都没有解析解，因此发展出各种近似 方法。伽辽金加权余量法与罩兹( r i t z ) 变分法是早期求解电磁场控制方程( 2 2 1 ) 0 0 两种 最常用的近似方法，也是建立有限元方程最重要的方法。 6 硕十学位论文 基于三棱柱嘲格的三维矢量有限元方法研究 假定面是控s f j y y 币1 ( 2 2 1 ) 的近似解。将西代入式( 2 2 1 ) ，由于与真实解之间存在 误差，因此会得到一个非零余量： ，= 三面一g ( 2 2 6 ) 伽辽金加权余量法通过对上述余量求加权积分： 5j a w ，r d q ( 2 2 7 ) 来寻求式( 2 2 1 ) 的近似解面，使得参数，在所有点上取得最小的值。上式中的m 为加权函数。近似解西可以表示为下面的展开函数形式： 面= q z ( 2 2 8 ) 这里，称作插值函数，一般为线性或高次函数。权函数的选取有多种方法，如点 配置法、最小二乘法等，但在加辽金加权余量法中，权函数w 与近似解展开所用的 函数，相同，通常这样可得到最精确的解。 里兹( r i t z ) 变分法，是一种用变分表达式( 泛函) 来得到近似解的方法。这种方法将 边值问题用泛函表示，泛函的极小值对应于给定边界条件下的控制方程的解。( 2 2 1 ) 式的近似解西通过求下面的泛函对面的极小值来得到： f ( 面) = 兰( 三面，面) 一三( 面，厂) 一三( 厂，西) ( 2 2 9 ) 方程( 2 2 9 ) 的极值点对6 f = 0 对应于( 2 2 1 ) 的近似解。因此一旦确定了方程( 2 2 1 ) 对应的泛函，即可将面的展开式代入来求解。 一般来说，用伽辽金加权余量法或旱兹变分法导出的方程具有同样的形式。对上 述基于电场的矢量h e l m h o l t z 波动方程( 2 2 2 ) ，得到的矢量有限元方程如下： 工去( v e ( r ) ) ( v 州r ) ) 一碍。e ( r ) w ( r ) d v = - j k o z o j ( r ) 州( 2 2 1 0 ) 上述方程是对整个问题区域进行建模得到的。如果所分析的问题过于复杂，则不 能使用简单的函数展开来近似整个区域内的解，就需要对整个问题区域进行离散，并 在每个子域使用伽辽金或里兹变分法来得到整个问题的解，这就是下面要介绍的有限 元方法。 7 2 有限元法的牡奉原理 硕士学位论文 2 2 3 有限元法的步骤 应用有限元法求解电磁场边值问题，一般包含下面的几个步骤： ( 1 ) 确立适当的微分控制方程及边界条件。 ( 2 ) 网格离散。 ( 3 ) 选择基函数和加权函数，运用伽辽金加权余量法或里兹变分法将控制方程离 散为线性方程组。 ( 4 ) 消去边界上的未知量并求解矩阵方程，得出所分析区域内的场分布。 ( 5 ) 后处理，计算出所需的参数。 使用有限元法，一旦确定了控制方程及计算求解区域，就需要进行网格离散。在 有限元分析步骤中，区域离散是非常重要的一步。因为网格离散的方式将会影响到计 算机内存需求、计算时间和数值结果的精确度。有限元法可以适应不同类型的网格， 对二维问题可以采用三角形、矩形及任意四边形等网格；而对三维问题可以采用四面 体、四棱锥、三棱柱、六面体等网格剖分。目前最为常用的，在二维问题中为三角形， 在三维问题中为四面体网格。正是离散网格的多样性使得有限元法在模拟任意形状边 界的问题时具有强大的功能。 图2 2 1有限元四面体单元示意图 在生成网格之后，就要根据选定的插值函数来建立有限元方程。对于切向矢量有 限元通常采用w h i t n e y 基函数，因为其由w h i t n e y 提出而得名。w h i t n e y 基函数形式 如下： w ，= ( “v 0 。 显然， 删母锯曷 2 舢 为了简化标记，也为了与在全局坐标系统，即直角坐标系下的坐标z 相对应，本 文把上述一维线性单元的局域坐标记为f ，而f 就是独立坐标六，1 - ( = 磊也就是原 来所讨论的非独立坐标。也即 3 基于三棱柱嘲格的三维矢量有限厄方法 硕士学位论文 其中， pz z l 2 。 甩 h = z 4 一z 1 2 、参数坐标舌，色和磊x 平面上的三角形面积坐标 ( 3 2 5 ) ( 3 2 6 ) 因为直角坐标系下三棱柱单元的上、下底面平行于x o y 平面：而在局域坐标系下， 标准三棱柱单元的上、下底面平行于卣c 彘平面，所以，本文通过一种较为简单的途 伞三 2 m = 札( t 乃- x ，y 2 ) + x ( y 2 - y 3 ) + y x 3 - x 2 ) 8 = 三任兰蔓l = 三c 而弘+ 五蜴+ 而乃一而乃一而躬一而弘， c 3 2 8 ， 1 6 3 2 图3 2 4 二维自然坐标 为便于标识，定义如下记号： q = 而弘一y 2 x a ；b l = y 2 一y 3 ；q = x 3 一x 2 口2 = x 3 y i y 3 x l ；6 2 = y 3 一y l ；c 2 = 而一x 3 ( 3 2 9 ) g 3 = x l y 2 一y l x 2 ；b 3 = m y 2 ；c 3 = x 2 一耳 硕士学位论文 基于三棱柱网格的三维矢量有限兀方法研究 于是，a 1 司以写成 l = 去( 口l + b l x + c j y ) 因此， 磊= 会= 去( 口1 + 6 l x + q y ) ( 3 2 1 0 ) 同样可得到 ， 乞= 叁= 西1 ( 6 2 x + c 2 y ) ( 3 2 1 1 ) 以及 一 磊= 叁= 去( q + 6 3 x + 巳y ) ( 3 2 1 2 ) 式中，：表示由点尸和结点3 以及结点1 组成的三角形面积：同样，表示由 点p 和结点1 以及结点2 组成的三角形面积。 因为点尸的选择完全定义了石，参和磊，所以，本文反过来，可以用它们来定 义三角形内的一个点。因此，( 当，岛，磊) 也被称为面积坐标。 妻 = 互 差差兰 三) i磊+ 岛+ 磊= 1 点五+ 受而+ 乞毛= 工 【卣乃+ 参+ 六儿= y 所以，用上面定义的面积坐标，容易建立一般三角形单元的结点插值函数。 可见，局域坐标系下的结点坐标具有如下特点： 平行于端点为2 、3 的线段的点的参相等，其他亦然。 处于端点为2 、3 的线段上的点的卣= 0 。 顶点1 的坐标为l ，即毒= 0 。 1 7 3 基于三棱柱网格的三维矢量有限冗方法 硕士学位论文 表3 2 3 局域坐标系下三角形各顶点的面积坐标 莓迫 磊色乞 d l loo 见 0 lo qool 即日( 1 ，0 ，o ) ，d 2 ( o ，1 ，o ) ，b ( o ，0 ，1 ) 经过上面的推导，容易发现，结点基函数具有以下性质： 3 、本部分结论 1 8 磊( ) 母$ 1 = j l j 图3 2 5 三棱柱单元局部坐标系 彘 硕十学位论文 基于三棱柱网格的三维矢量有限元方法研究 表3 2 4 三棱梓单元各顶点的整域坐标及与之相对应的局域坐标 坐标 整域坐标局域坐标 葡八 xyz 量毒2 l五j ，l z l l00 2 屯儿 z l olo 3 j c 3乃 z 1 o0 o 4 mz 4 10l 5 x 2y 2z 4 oll 6 x 3y ，z 4 ool 建立整域坐标和局域坐标之间的映射关系： 整域坐标系下，三棱柱单元的母线平行于z 轴，上、下底面平行于x o y 平面，而 且底面和顶面的三角形的棱边不必平行于x 轴和y 轴。 局域坐标系下，卣，参和f 是独立坐标，与是非独立坐标；规则三棱柱单元的 母线平行于f ，上、下底面平行于缶c 彘平面；在底面三角形上f = 0 ，在顶面三角形 上f = 1 。 磊，受和岛是二维平面三角形的面积坐标，f 是一维线性单元的位置坐标，并 且有 点：i 丢( q + 6 l x + q y ) 向2 西+ 岛x + q y j 参= 去( 吒脚坳) ( 3 2 “) 岛= 去( 巳+ 岛石+ 巳y ) 尸一z z l ( - = o 其中，单元插值函数系数口，、岛、q ，i = l ，2 ，3 的计算表达式见( 3 2 9 ) ； 单元体的母线长度h = z 。一z 。： 1 9 3 基于三棱柱嘲格的二维矢量自限元方法 硕 ：学位论文 单元上、下底面三角形的面积8 的计算表达式见( 3 2 8 ) 。 这些参数坐标之间的关系及各自的数值范围是 缶+ 岛+ 岛= l ( 3 2 1 5 ) 且 卣，参，岛，f 【0 , 1 】 ( 3 2 1 6 ) 3 2 1 3 在局部参数坐标系下进一步讨论单元的几何表示 定义是指向单元内部的位置矢量。矢量是从直角坐标系下的实际单元到参数 坐标系下由一组独立参数坐标定义的标准单元变换过程中得到的参数映射。这组独立 的参数坐标本文在前面定义局部参数坐标系的时候已经提及，它们是( 卣，色，f ) 。 1 、单元棱边矢量 为了定义三棱柱单元的棱边矢量基函数，首先定义单元棱边矢量： 1 1 ：o r1 2 ：旦1 4 ：o r( 3 2 1 7 ) a 专、a 专、8 雅可比行列式 j = 0 7 1 7 ) 1 。 ( 3 2 1 8 ) 其中， ，k 的取值为 1 ，2 ，4 ， 2 ，4 ，1 或者 4 ，1 ，2 。 当三棱柱是直线型三棱柱时，g l , e 2 , 粤4 是常量矢量。棱边矢量粤与参数坐标的梯 度存在双正交性的关系： i j 飞毛l = 6 q f 1 ，i = j f ( 3 2 1 9 ) 1 0 ， i _ 由上述双j 下交性关系，对三棱柱单元，建立下面的关系： v 专= 了fx !k(3220) 也就是说， v 卣= 等v 磊= 等v f = 了! i x1 2 ( 3 2 2 1 ) 硕l 学位论文 堆于三棱柱刚格的三维矢量有限元方法研究 凼此，对独立参数坐标来说，同样司以建立这样的关系： v 乞x v 磊= 专1 7 ( 3 2 2 2 ) 显然，棱边矢量的方向，指向一个梯度矢量v 乞和另一个梯度矢量v 磊的叉乘所 得的矢量所指的方向。 三棱柱单元中，每条棱边都对应一个棱边矢量，于是每条棱边的棱边矢量存在如 下关系： i 2 = 1 5 = i 1 ：1暑u=i21 1 2 ， 1 3 26 = 1 2 1 1 7 = 1 8 = i 9 = 1 4 最后，给定了三个一组的非共面参数坐标变量卣，色和f ，有 v 缶( v 乞v f ) = 专( 3 2 2 4 ) 该雅可比行列式严格为正，这也是之前建立局部参数坐标系的原则之一。 2 、单元梯度矢量 定义参数坐标变量的梯度如下： v 参一每 2 筋， 也就是说，定义坐标变量梯度的大小为一1 ，而梯度矢量的方向则与一个相应的 r l 单位矢量蓐所指向的方向相反。单位矢量蓐指向单元的外部。 一组参数坐标中，独立坐标的梯度矢量和非独立坐标的梯度矢量存在下面的关 系： v 岛= 一v 卣一v 岛 ( 3 2 2 6 ) 参数坐标的梯度矢量 2 l 3 基于三竺壁塑整塑三笙叁垦鱼：堡：i 查鲨堡主兰竺兰l 梯度矢量的叉乘 v 卣= 薯；+ - 粥f l y - = 去( 岛；+ q 多) v 受= 警；+ 等多= 去( 6 2 ；均纠 2 刀， v 磊= 警王学多= 去( 6 3 ；+ 岛多) v 彳：丝o z 三= 去三 第一类通式，是v 专xv g ，共甲，卢l ，2 ，3 ： 根据上面求得的表达式v 喜= 去；+ 否e ly a ，v f = 去三，所以， v 毒v f = 去( q 三一岛多) 即 v 岛v f ；去( c ：；一6 2 多) v 磊v f = 去( 巳；一6 3 多) v 卣v f = 去( q ；一6 i 多) ( 3 2 2 8 ) ( 3 2 2 9 ) 第二类通式，是v 毒v 乞，其中，扛2 时，= 3 ；或者待3 时，j = l ；又或者 i = l 时，_ ，= 2 ： 即 v 触乞2 卉( 盱q 屯) 三 ( 3 2 3 0 ) v 受v 磊2 研1 ( 6 2 2 6 3 z v 岛v 磊2 研1 6 3 6 1 z 3 2 3 1 v 铋2 卉n c 2 - c 1 6 2 乒 硕上学位论文 基于三棱柱嗍格的三维矢量有限元方法研究 ( d 梯度矢量的点乘 第一类通式，是v 专v 乞，其中，i = 1 ，2 ，3 = 1 ，2 ，3 ： v 专v 乞2 土2 a 生2 a + 丢丢2 卉( 勿q + q 刁) 2 抛， 第二类通式，是v f v f ： v ( o r f = 吉 第三类通式，是v 毒v f ，其中，i = 1 , 2 ，3 ： 3 、单元高度矢量 v 毒v f = 0 虿= 曩五= 一曩2 v 毒 ( 3 2 3 3 ) ( 3 2 3 4 ) ( 3 2 3 5 ) 4 、本部分小结 下图所示为三角形单元和三棱柱单元的棱边矢量、高度矢量和梯度矢量。 | 沙一 i 广 7 、v 号3 v 亏。 f j 1 一 图3 2 6 三角形单元的棱边矢量、高度矢量和梯度矢量 3 幕于三棱柱| 旬9 格的三维矢量自限几方法 硕十学位论文 图3 2 7 三棱柱单元的棱边矢量、高度矢量和梯度矢量 3 2 2 单元矢量基函数的选择 如果单元每条棱边被赋予一个常切向分量，那么该单元内的矢量场可以展开为： e 。= n 露 ( 3 2 3 6 ) 其中，群表示沿着第f 个棱边的切向场。吖是三棱柱单元的矢量插值函数或基 函数。 首先，对三棱柱单元的各个结点和棱边进行局部编号，然后，引出相应棱边的矢 量基函数。对每个三棱柱单元，都有9 个矢量基，图示如下： 硕j 二学位论文基于三棱柱网格的三维矢量有限元方法研究 下面给出三棱柱的棱边矢量基函数的公式 n 。= ( 1 一f ) ( 受v 乞一磊v 受) n ：= ( 1 一f ) ( 磊v 卣一点v 磊) n 3 = ( 1 一f ) ( 卣v 受一磊v 卣) n 。= f ( 磊v 磊一乞v 乞) n ，= f ( 岛v 缶一卣v 磊) ( 3 2 3 7 ) n 。= f ( 磊v 受一色v 卣) n ，= 考翼 n s = 毛羹乞 n ，= 芎汴 其中， 三棱柱单元上的上、下底面三角形单元的矢量基函数选取的是w h i t n e y 1 型基函 数，即如= 氧v 乞一专：v 4 。石、岛和夤是二维三角形单元的面积坐标，也就是结点 基函数。 f 位于上、下底面三角形之间，是一阶拉格朗日插值多项式。当位于三棱柱单元 的底面时，f = 0 ，当位于三棱柱单元的顶面时，f = 1 。 所以，基函数n 、n ：和n ，位于三棱柱单元的底面上，它们是由参数坐标( 1 一f ) 乘以二维三角形单元的棱边矢量基函数w 得来的； 基函数n 。、n ，和n 。位于三棱柱单元的顶面上，它们是由参数坐标彳乘以二维 3 基于三棱柱 6 9 格的三维矢量有限元方法 硕上学位论文 三角形单元的棱边矢量基函数形得来的： 基函数n ，、n 。和n ，与三棱柱上、下底面相垂直，与之相关的未知量被称作体 未知量。 把三棱柱单元的9 个棱边矢量基函数分情况表示如下： 底面矢量基 顶面矢量基 体矢量基 n = ( 1 一f ) w 2 ， n z = ( 1 - f ) w 3 - ( 3 2 3 8 ) n ，= ( 1 一f ) w l ： n 4 = f w ：3 n 5 = f w 刍 n 。= f w l 2 n 1 = 考 n 8 = 参v f n ，= 专孓专 ( 3 2 3 9 ) ( 3 2 4 0 ) 其中，三角形棱边矢量基函数是 w & = 氧v 彘一专：v 氧 ( 3 2 4 1 ) 具体有 w 凸= 考零考3 一专孓专2 w j i = 磊v 点一点v 磊 ( 3 2 4 2 ) w l ：= 当v 彘一磊v 磊 这样建立的棱边矢量基函数满足矢量基的基本特性，即沿相关棱边有常切向分 量，而沿其他棱边则无切向分量。 3 2 3 单元s 矩阵和t 矩阵的生成 当将上面建立的矢量基函数用于表示矢量波动方程有限元解中的矢量场或势时， 有必要考虑与之相关的单元矩阵的计算。从矢量波动方程离散而产生的单元矩阵包含 下面两种形式的积分： 鬈2 罂r ( v 刚( 乳m 。p y ( 3 2 4 3 ) f ：= 髓i | n ：n ；d v 以三棱柱单元为剖分单元，必须对这两个积分进行数值计算。 硕上学位论文基于三棱柱网格的三维矢量有限元方法研究 3 2 3 1 三棱柱单元各个棱边矢量基函数的旋度 1 、底面矢量基的旋度 v n t = 乞v 毒v 0 一f ) 一氧v 气v ( x f ) + 2 ( 1 一f ) v 氧v 气( 3 2 4 4 ) 因为v ( 1 一f ) = - v f ，所以， v n i = 氧v 考：v f 一乞v 氧v f + 2 ( 1 一f ) v 氧v 气 ( 3 2 4 5 ) 代入梯度矢量叉乘公式( 3 2 2 8 ) 式和( 3 2 2 9 ) 式， v n i = 去 ( 话喵z ) ；一( 盹乞) 多 + 卉( ”砌m 乒 ( 3 2 4 6 ) 简记为 v n ，= ( 厶氧一乞) 三+ ( 蜀乞一g t 善1 1 ) 多+ b ( 1 - f ) 三 ( 3 2 4 7 ) 1w h e n i = 1 ， = 2 ，f 2 = 3 ； 其中， w h p 疗i = 2 ，毛= 3 ，之= l ； 【w h e n i = 3 ， = l ，f 2 = 2 ； 2 、顶面矢量基的旋度 很明显，对三棱柱单元来说，它的底面矢量基和顶面矢量基之间的差别只是f 因 子：三棱柱单元的底面矢量基，是由参数坐标( 1 一f ) 乘以二维三角形单元的棱边矢量 基函数而得：