（电工理论与新技术专业论文）三维静电场高精度边界元方法的研究.pdf

声明 = 牟：人郑霞卢明：此处所提交的硕+ 学位论文二维静电场高精度边界元方法的研 究，是本人在华北电力大学攻读硕士学位期间，在导师指导下进行的研究- t 作和取得 的研究成果。据奉人所知，除了文中特别加以标注和致谢之处外，论文中小包含其他人 已经发表或撰写过的研究成果，也小包含为获得华北电力大学或其他教育机构的学位或 证1 亏而使用过的材料。与我一同t 作的同志对奉研究所做的任何贡献均已在论文中作了 i i l j 确的说明并表示了谢意。 学位论文作者锯名：翌查暨 日期：巫三堡 关于学位论文使用授权的说明 本人完全了解华北电力大学有关保留、使用学位论文的规定，即：学校有权保管、 并向有关部门送交学位论文的原件与复印件：学校可以采用影印、缩印或其它复制手 段复制并保存学位论文：学校可允许学位论文被查阅或借阅；学校可以学术交流为 目的。复制赠送和交换学位论文；同意学校可以用小同方式在不同媒体上发表、传播学 位论文的全部或部分内容。 ( 涉密的学位论文在解密后遵守此规定) 作者签名： 等 华北电力人学硕十学化论文 第一章引言 1 1电磁场数值计算方法概况 电磁场理论的研究具有悠久的历史，而计算电磁场使整个电磁场理论发生了苹 命性的变化，计算电磁场对电磁场理论发展的影响决小仅仅是提供一个计算t 具而 已，日前，计算电磁场已成为对复杂体系的电磁规律、电磁性质进行研究的手段， 为电磁场理论的深入研究开辟了新的途径，并极大地推动了电磁场t 程的发展。 常用的计算电磁场问题的方法主要有两大类，其中每一类又包含若于种方法， 第一类是解析法：第二类是数值法。数值法的出现，使电磁场问题的分析研究从解 析的经典方法进入到离散系统的数值分析方法，从而使许多解析法很难解决的复杂 的电磁场问题，有可能通过电磁场的计算机辅助分析获得很高精度的数值解，同时 极大地促进各种电磁场数值计算方法的发展”。 已经发展起来的数值计算方法可以分为微分方程法和积分方程法两大类。微分 方程法主要包括有限差分法( f d m ) 和有限元法( f e m ) ；积分方程法可再分为模拟 电荷法( c s m ) 和边界元法( b e m ) “1 。 在电磁场数值计算方法中，有限差分法是廊用最甲的一种方法，它是以差分原 理为基础的一种数值计算方法。2 0 世纪5 0 年代以来，有限差分法以其数学概念清楚、 形成系数矩阵f 分方便等特点，在电磁场数值分析领域内得到了广泛的应用。但有 限差分法的规则嘲格一i 能满意地模拟几何形状复杂的问题，而在电t 设备中的电磁 场却往往是以包含复杂的几何形状和小同材料的物理参数为特征，因此有限差分法 在电磁场分析中的应用逐渐被有限元法替代”。 有限元法是以变分原理和剖分插值为基础的一种数值计算方法。在力学领域。 有限元的思想早在2 0 世纪4 0 年代已经提出，在2 0 世纪5 0 年代开始用于飞机设计。但 是这个方法开创性的工作，公认是r w c l o u g h 在1 9 6 0 年发表的著作中奠定的”。此 后，该方法得到发展并广泛地用于结构分析，流体力学、热传递等物理和t 程问题 中。 2 0 世纪6 0 年代末至7 0 年代初，有限元法被移植到电磁场丁程领域的设汁汁算 中。1 9 6 5 年，w i n s l o w 河先将有限元法应用于电气- t 程问题”1 ，其后，1 9 6 9 年s i l v e s t e r 将有限元法推广应用于时谐电磁场问题”。7 0 年代a n d e r s o n 对变压器漏磁场的研究， o k u d a 等人对气轮发电机端部磁场的研究，8 0 年代n a k a t a 等人对电磁材料特性的数 值模拟和实验研究都是富有开创性的成果。目前，有限元法已经在电磁场数值汁算 方法中占据主要地位 模拟电荷法是一种等效场源法。最早由h s t e i n b i g l e r 提出。模拟电荷法原理 简单明了，在具备场的基奉概念的基础上，解出离散的模拟电荷值后，场量即可由 华北电力人学硕士学化论文 解析式算得，有较理想的场量计算精度，冈此成为静电场数值计算的主要方法之一 ”。但当电极形状或场域边界很复杂时，模拟电荷的数量往往取得很大，模拟电荷 法也必须采用系统的处理手段使问题能够易于在汁算机上求解，所以在一定程度上 限制了模拟电荷法的应用。 边界元法是在经典的边界积分方程法的基础上发展起来的一种数值计算方法。 它可以分成两种基奉类型，即间接法和直接法”1 。间接边界元法是用物理意义小一 定很明确的变量来表示化成的公式。这种方法曾用于求解由拉普拉斯方程或亥姆霍 茨方程所控制的弹性力学问题或其它势问题”1 。直接边界元法是上世纪8 0 年代提出 来的，它是以格林恒等式作为出发点，变量具有明确的物理意义。现在已优先府用 于t 程科学中，= 牟：论文仅讨论研究这种方法，这一方法在我国也引起了许多学者的 晕视，如冯康呻1 ，牡庆华1 、何广乾、劳洁生、曲圣年等。日前，边界元法已经成 为电磁场数值计算方法的主要方法之一。 现在，还出现了把上述几种电磁场计算方法相结合起来的新方法，即耦合法”“。 目前，用来处理电磁场问题的有有限元与边界元、有限元与模拟电荷法、有限元与 积分方程法互相组合的方法。因为耦合法可以取长补短，所以越来越得到大家的关 注 综上所述，当前主要的电磁场数值计算方法有有限元法、边界元法和耦合法。 同时也出现了多种成熟的电磁场数值分析商业软件，如a n s y s 、e l e c t r o 、c o u l o m b 、 m a f i a 等。 1 2 三维静电场边界元法精度研究的必要性 电场计算是高压电气设备设计中的霞要内容之一，一个合理的高压电气设备设 计，小仅要求最大电场强度低于相应介质的介电强度，而且要求电场分布尽可能均 匀。 在高压电气设备中导体的附近存在较高的电场强度，冈而对其绝缘介质提出了 较高的要求，例如导体引线的绝缘套管、导体相间和对地绝缘等。正确合理地选择 绝缘介质的材料、尺寸和形状，决定于电场强度计算是台正确。高压电气设备中相 问距离和对地距离的确定，也决定于最大电场强度的汁算”。在电力系统中，合成 绝缘子的电场计算更具有实际意义。合成绝缘子用于支撑、隔离或包容高压带电导 体。由于在线路运行中会出现各种各样的问题而导致绝缘子性能下降无法正常运 行，所以线路绝缘子故障是线路发生故障的主要表现形式之一，冈此要研究和改进 合成绝缘子的性能就很有必要对合成绝缘子进行电场分析计算。 对于合成绝缘子电场分析模型，是个开域场的问题。而有限元法的主要缺点之 一就是不大适合求解无限边界场域边值问题，所以用边界元法求解三维静电场的开 域问题就显得i 一分重要。 2 - 华北电力人学硕+ 学位论文 由于三维静电场边界元法通常采用近似的数值积分法，要怨获得较为精确的数 值解就需要我们去研究和总结提高计算精度的算法。 1 3 本文研究的主要内容 边界元法( b e m ) 是适合于求解电磁场开域问题的一种方法。这种方法是把所 研究问题的微分方程变成积分方程，然后将区域的边界划分为有限个单元，也就是 把边界积分方程离散化，得到只含有边界上的结点未知量的方程组，然后进行数值 积分求解。边界元法在单元上所考虑的函数可以按巧i 同的形式变化，这一做法与有 限元法大致相同。冈此，在三维静电场边界元法分析中一个霞要的课题就是选择合 适的单元和单元内合适的插值函数来提高边界元法的计算精度。 日前已有的三维静电场边界元法分析中，均采用的是线性等参元有限元单元 ”1 。如果所求边值问题需要更高计算精度，就需要具有更高阶次的插值函数 1 7 1 ， 也就需要高阶等参元有限元单元，所以本论文编制了三维静电场边界元法的二阶等 参元有限元单元的计算软件。 奉论文编制的计算软件，是以有限元分析软件a n s y s 进行前处理即进行建模和 嘲格剀分，采用v i s u a lf o r t r a n 6 5 编写的。 此外，本论文还主要进行了以下几个方面的下作： 1 ) 研究了三维静电场边界元法二阶单元的积分精度问题。分析影响该积分精度的 相关冈素，即单元形状、大小、个数和高斯点个数等对积分精度的影响。分析 如何改善积分精度。分析如何解决奇异积分的高斯积分误差。 2 ) 介绍基于球坐标三维静电场曲边四边形边界元方法和基于圆环坐标的三维静电 场曲边四边形边界元方法。 3 ) 利用三维静电场边界元法线性单元、二阶单元程序和基于球坐标及基于圆环坐 标的曲边四边形边界元程序分别对不同的电极模型进行计算，对比分析这几种 边界元算法对计算精度的影响。 4 ) 从理论和实践上充分比较三维静电场边界元法的四种提高计算精度方法的优缺 点，总结规律，为以后的三维静电场边界元法计算提供参考依据以便后人选择 使用。 华北电力人学硕十学化论文 第二章边界元法基本理论 2 1引言 边界元法( b o u n d a r ye l e m e n tm e t h o d ) 是从2 0 世纪7 0 年代发展起来的一种数 值计算方法。该方法的t 程应用起始于同体力学，习前在流体力学、热力学、电磁 学、土木t 程等领域也得到广泛的应用，并已从线性、静态问题延拓到非线性、时 变问题的研究范畴。 边界元法是把边值问题的微分方程变成边界积分方程，然后将区域的边界划分 为有限个单元，也就是把边界积分方程离散化，得到只含有边界上的结点未知量的 方程组，然后进行数值求解。有限元法和边界元法都是离散化数值解法，它们有各 自的优缺点和局限性。 边界元法有如下的主要特点： 1 ) 降低问题求解的空间维数。边界元方法将给定空间区域的边值问题通过包 围区域边界面上的边界积分方程来表示，从而降低了问题求解的空间维数， 即三维问题可利用边界表面积分降为二维问题；而二维问题则利用边界的 线积分降为一维问题。冈此，有限元离散仅对廊于二维曲面单元或一维曲 线单元，大大简化了问题的前处理过程。 2 ) 方程组的阶数降低，输入数据量减少。由于待求的结果仅限于边界结点， 所以大幅度降低了离散方程组的阶数，节省了计算机的内存，同样，简化 了问题的前处理过程。 3 ) 计算精度高。除了对狭长形状的解域外，边界元法的求解精度一般高于有 限元法。边界元法直接求解的是边界广义场源的分布。根据4 i 同的问题， 广义场源可以是位势、场源或等效场源。场域中任一点的场量可通过线性 叠加各离散的广义场源的作用而求得，卅i 需要再用微分运算求得。另外， 由于只对边界离散，系统误差仅仅来源于边界，所以边界元法有较高的计 算精度。 4 ) 适于求解开域问题。有限元法的主要缺点之一是不大适合求解无限边界场 域边值问题，而边界元法可对有限场域和无限场域的有限边界进行离散化 处理并求解。 5 ) 不易处理多种媒质共存的问题。由于边界元法在场域中利用了自由空问的 格林函数，因而4 i 适合求解含有非均匀媒质的边值问题。日前出现的耦合 法能扬长避短解决许多复杂的工程问题。 2 2 边界元法基本理论 华北电力人学硕十学化论义 2 2 1 标量格林定理 设某一闭合曲面s 所罔区域为矿。若有两个标量函数p 和。它们在v 域内及s 面上分别有连续的一阶和二阶偏导数。根据矢量恒等式有 v ( 妒妒) = 卯2 矿+ v 妒v p ( 2 - 1 ) 根据散度定理有 j m 妒2 y + 咖( v y ) 】d y = 辑( 椤) - d s ( 2 2 ) vj 式( 2 - 2 ) 称为格林第一恒等式。将妒与交换可得 肌汐2 妒+ ( v ) 一( v 例d y = 彤( 咿力d s ( 2 3 ) r5 将式( 2 - 2 ) 与式( 2 - 3 ) 相减可得 j j 刃2 一咿2 o d v = 西( 妒一咿妒) - 拶 ( 2 4 ) v5 武( 2 - 4 ) 称为格林第二恒等式，又称为格林定理。 2 2 2 静电场的边界积分方程 获得边界积分方程有三条途径：叠加原理、等效原理即格林定理和互易性。奉 文边界元积分方程是通过格林定理推导出来的。 等效原理指区域v 内的源在k 中产生的效应或k 内的源在矿内产生的效麻，可 用形和矿两区域的闭合分界面上的等效源的面积分等效替代，被替代区域的内部场 量为零0 7 。 如图2 - 1 所示，设有电荷体密度p 分布在区域v 内，p 1 分布在区域k 内， 若只需计算体积v 以内区域中任意点 p 点的电位，应是场源p 和岛对v 内备 点产生的影响。先考虑p 对v 内的影 响。 静电场的基奉方程为 一 k 图2 - 1p 和n 对矿内并点产生的效戌 一田2 妒= p ( 2 - 5 ) 求解区域为矿，求解区域边界为s ，矿为v 内的电位，p 为矿内分布的电荷体 密度。假设求解区域内充满同一均匀介质。 在三维情况下，取自由空间的格林函数为 华北电力人学硕十学佗论文 炉习匆 他。6 式( 2 - 6 ) 中，和r 分别表示场点和源点到坐标原点之间的距离。由于场点和源 点的距离有可能为零，即有可能出现奇点，所以在整个区域内满足 v 2 缈= 一声( ，一，) ( 2 - 7 ) 式( 2 - 7 ) 中j 为狄拉克函数，应用格林定理可得 l 南c 一争删，卜亭c 南鲁一伊磊0 羽1 郴协s ， 根据狄拉克函数的性质，上武麻为 加，= l 詈南+ g c 豳鲁一号未向搬 协。， 另外，由叠加原理可得 舯f 詈南+ i i 譬尚( 2 - 1 0 ) 比较式( 2 9 ) 和( 2 1 0 ) 可知，式( 2 - 9 ) 中的面积分这项所求等效于形中p l 在v 内产生的效应，也就是可用包同v 的边界闭合面上的等效面电荷( 密度为0 q 0 1 0 n ) 和等效电偶极子层( 密度为印) 替代区域v 之外客观存在的场源。这是奉文计算软件 中边界元算法的根= 牟= 依据。 在求解区域边界面为光滑曲面的情况下，若场点落在边界面上，冈场域外c a 位 为零，则有 扣2 詈孝和相等一考嘉南，嬲( 2 - 1 1 ) 式( 2 - 1 1 ) 左边含有边界上的电位，右边含有边界上的电位和边界上电位的法 向导数。这就构成了边界电位及其法向导数的积分方程，简称边界积分方程。 如图2 - 2 所示，在求解区域边界上有 丽0 网1 = 杀c = c v 巳= 去印q = 一古c 。s 护 cz 一z ， 式( 2 1 2 ) 中，p 是和边界咖外法线方向气之间的夹角。因此边界积分方程又可 表述为 扣，= f 篆+ g c 去等+ 丢警，嬲 c z 一 华北电力火学硕士学位论文 图2 2 边霁萄上潦点与场点豹关系 以上是针对三维静电场i 司题推导的边界积分方程。对于二维静电场问题，令格 林函数为 妒= 如去 浯 通过与三维静电场炎似的推导，可得边界积分方程为 圭舯) = 一- 景8 i n r d s f ( 去n 嗉一丢警灯 c z 一- s ， 在二维求解区域边界上有 百8 1 n r = v ( 1 1 1 r ) 巳= 一i 1 巳= 去c o s 口 ( 2 一1 6 ) 式( 2 - 1 6 ) 中，口是- - e 。, 和边界面外法线方向巳之间的夹角。因此边界积分方程又可 表述为 三肿) = 一f 去l n r d s f ( 去l n r 警一瓦t p c o f s g 珂 ( 2 一1 7 ) 式( 2 1 3 ) 和式( 2 1 7 ) 即为用加权余量法进行离散化的边界积分方程。 2 2 3 加权余量法 2 2 3 1加权余量概念 假定有一边值问题满足方程 l u = g ( 2 - 1 8 ) 式中，“为未知函数；上是算子，表示对“的一种运算，既可以是微分运算，也可为 积分运算；g 为已知函数。 为了求出“，设有一组线性无关的函数，“：，以，取其前h 项的线性组合作 为“的近似解“。若当n 一时，有“_ “，则称，“2 ，以，为基函数序列，l i l t 为基 函数。“的近似解表示为 华北电力人学硕十学位论义 五= o 叶 ( 2 1 9 ) j - i 武中q ，巳，巳为待定系数a 将近似解代入方程得余量为 r = l u g = 艺巳砒- g ( 2 2 0 ) j - l 如果余量为零，说明已经满足方程，即一u 是方程的准确解。但一般情况下无法 做到余量为零，冈此只能放松约束强制余量的加权积分为零，即 l w f 肘q = 0( i = i ，2 ， ) ( 2 2 1 ) 式中为权函数，w 2 ，w k ，为权函数序列，权函数之间要求线性无关。权 函数的不同选择导致不同的近似方法。本文边界元法采用的是伽辽金法。 2 2 3 2伽辽金法 选取的权函数序列与基函数序州相同，即w t = 珥，此处的坼与式( 2 1 9 ) 中“，相 同t 将m 代入式( 2 - 2 1 ) 可得 l 坼r d f l = 0 ( f = 1 ，2 ，”卅) ( 2 2 2 ) 将式( 2 - 2 0 ) 代入式( 2 - 2 2 ) 可得 l 坼三杰勺哟d q = l ，u , g e n ( f = 1 ，2 ，) ( 2 2 3 ) 式( 2 - 2 3 ) 即伽辽金法的方程。其中的n 表示箅子的定义域，它可以是三维空 间、空间曲面，平面及平面的的边界线等。n 为离散结点的总数。a 其实是结点处 的函数值，而玑作为基函数，它是场点坐标的函数。 由于c 是结点的函数值而不是坐标的函数，所以式( 2 - 2 3 ) 可改写成为 喜勺叫q = “枷 上式可用矩阵形式表示为k c = 6 。这样，通过加权余量法的伽辽金格式我们把 算子方桓转化成了代数方程组。当基函数被确定后，印系数矩阵k 和右端项矩阵6 被 确定，我们就可计算出定义域中各结点的函数值 2 2 4 插值函数与形函数 有限元法和边界元法的首要问题就是确定特定的单元和单元内合适的插值函 数。这个问题进行得好坏将直接影响到计算结果的精度，冈此，它是本论文研究的 出发点。 虽然许多类型的函数都可作为场变量的插值函数，但拉格朗甘插值多项式是廊 8 华北电力人学硕十学化论文 用最为广泛的。主要原因有两个；首先。在建立备类单元的方程和数宁计算时，多 项式的数宁处理比较容易，尤其是微分和积分的运算比较容易；其次，任意阶次的 多项式可以近似地表示真实解。 在实际虑用中，对于一个具体的单元，用来表示场变量的插值多项武的阶次， 一般不能任意确定。多项式的阶次数取决于该单元结点变量的数日( 称为自由度) 。 换言之，多项武中系数的数日廊等于单元的自由度数。通过单元结点的函数值，可 以联立方程组求得多项式的系数。这样，插值函数就可表示为结点的形状函数和结 点函数值的乘积之和。 形状函数为场点坐标的函数，它的具体表示式与单元的几何形状和单元内的结 点数有关，故又称为单元形状函数。 2 2 5 整体坐标与局部坐标 在二维场域中备点的位置，既可采用直角坐标，又可采用极坐标，这两种坐标 可用于表示场域中所有点的位置，所以通常称它们为整体坐标。相对于整体坐标而 言，局部坐标只在每个单元中有定义，即它只用于描述每个单元中备点的位置，在 单元之外没有意义。自然坐标是一种局部坐标，它是用一组都小超过l 的无量纲数 来规定单元中点的位置，这种坐标系大大简化了计算，对平面单元和曲面单元同样 适用。 2 2 5 1 三角形单元的自然坐标 如图2 - 3 所示，设单元的面积为，单元内任一点p 与二个结点的连线将单元 分成三部分，它们的面积分别记作，。，定义p 点的自然坐标，又称面积坐标， 为 y j j 超 0 图2 3 三角形单元的面积坐标 厶= 会 = 等 。= 会 ( 2 2 5 ) 由于厶+ 厶+ o = l ，所以三个面积 坐标中只有两个是独立的，这样，三角 形单元内任一点的位置，只用厶，k 中任意两个表示即可。从面积坐标和整 体坐标的关系可知，面积坐标厶，k 分别 是三结点f ，j ，m 的形状函数如果指定厶= 0 的直线为f 轴，厶= 0 的直线为叩轴，如 图2 - 4 所示。则通过以上坐标变换，整体坐标中的任意三角形单元都可转化为自然 华北电力人学硕十学位论文 坐标中的直角边边长为1 的标准三角形单元。 圈2 4 三角形单元的0 一口坐标系 2 2 5 2 四边形单元的自然坐标 与三角形单元一样，除用整体坐标外，也常用局部坐标。通过坐标变换，可将 整体坐标中的任意四边形单元转化为自然坐标中的以原点为中心的边长为2 2 的标 准矩形单元，如图2 - 5 所示。同样，由自然坐标和整体坐标的关系可得到结点的形 状函数。 i 4 ( - 1 i ： 可 1 、 0善 1 ( 1 ，- 1 )文1 ，1 ) ， l 玎 4 r it 、1可1 、 ， 80； i ( 1 1 ) 5 瓦l ，1 ) 圈2 5 四边形单元的善一口坐标系 2 2 6 等参数单元 由上述可知，在场变量中用到了形状函数，而在坐标的变换式中也用到了形状 函数。这两种函数的性质相同，但选择的出发点却有所不同。对于前者，形状函数 的阶数越商，截断误差越小；对于后者，应根据场域边界形状的不同，选用相麻阶 数的形状函数，以使场域拟合得更好。 在实际使用的单元中，我们对场变量的形状函数和坐标变换式中的形状函数采 用了等同的形式，那么，这个实际单元就称为等参数单元即等参元。为了能以较少 的剖分单元获得较高的计算精度，使j | f | 分单元能与场域边界较理恕地拟合，同时为 华北电力大学硕十学位论文 了计算方便，本文采用的是二阶等参元插值边界元法。 2 3 本章小结 奉章介绍了边界元法的理论依据：标量格林定理、静电场的边界积分方程和加 权余量法。同时还介绍了场域离散中的插值函数与形状函数，引进了整体坐标与局 部坐标的概念以及等参数单元的定义，为第二章的三维静电场二阶插值边界元法的 研究打好了理论基础。 华北电力人学硕+ 学化论文 第三章三维静电场二阶插值边界元法 3 1 三维边界元算法 三维场对应的边界面单元为线性的三角形单元或四边形单元和二阶的三角形 单元或四边形单元。 3 1 1 单元类型 3 1 1 1 线性单元 在线性单元中，比较常见的就是具有3 个结点的三角形单元和4 个结点的四边 形单元。如图3 - l 所示。 2 图3 1 线性单元 2 在标准三角形单元上，3 结点的三角形单元的形状函数为 l = 1 一孝一叩，2 = 孝，3 = ，7 ( 3 - 1 ) 在标准矩形单元上，4 结点的四边形单元的形状函数为 3 1 1 2二阶单元 n i = 丢( 1 一# ) ( 1 一”) 2 = i 1 ( 1 + f ) ( 1 一玎) m = 丢( 1 + 善) ( 1 + 町) m = 丢( 1 一手) ( 1 + 野) 3 ( 3 - 2 ) 图3 - 2 二阶单元 对于高阶单元来说，本文主要研究的是二阶单元。在二阶单元中，比较常见的 - 1 2 华北电力人学硕十学化论文 是具有6 个结点的三角形单元和8 个结点的四边形单元。如图3 2 所示。 在标准三角形单元上，6 结点三角形单元的形状函数为 l = ( 1 一善一口) ( 1 2 4 2 r ) 乞= 善( 2 古一1 ) 3 = i ( 2 7 一1 ) 4 = 4 孝0 - # - r ) n s = 4 翻 n o = 4 1 ( i - 孝- r ) ( 3 - 3 ) 由于奉章重点介绍二阶插值单元，即三角形二阶插值单元的形状函数分别对 善，口求偏导得 嚣卅州峋等叫。 尝：o譬：4 喵却 a fa f 。 等却警= 一叩a a 。 丛：- 3 a 玎 当：4 ，7 d 行 擎：4 善 d 玎 + 4 孝+ 4 ，7 婴：o d 疗 一1丛：- 4 f 口玎 婴：4 4 f 一8 ，7 d 玎 在标准矩形单元上，8 结点四边形单元的形状函数为 j = ；( 1 一孝) ( 1 一，7 ) ( 一善一玎1 ) 3 = ；( 1 十善) ( 1 + 玎) ( 孝+ r - 1 ) m = ；( 1 管) ( 1 - q ) 7 = 兰( 1 一善2 ) ( 1 + ，7 ) 2 = 丢( 1 + f ) ( 1 _ ，7 ) ( f 一叩一1 ) = 去( 1 一手) ( 1 + ，7 ) + 刁一1 ) 6 = 圭( k ) ( 1 一矿) l = 吉( k ) ( 1 一町2 ) 同样，四边形二阶插值单元的形状函数分别对善，7 求偏导得 ( 3 - 4 ) ( 3 - 5 ) ( 3 - 6 ) 华北电力人学硕十学位论文 嚣= 扣堋训器= 扣训o - 玎， 等= 扣c 警= 扣训 等吲州，警= 扣相 等叫( 1 瓦a n s = l ( r 2 - 1 )d cd 0 鲁= 扣z 枞z 咱 警= 和z 坝，均 鲁= 抄叫 盟：三( 1 一尹) a t l 2 、 7 丛：三(2刁一,do+孝)at 4 、 77 婺：! ( 2 峭) ( 1 瑚4a 玎 、” 擎；一，7 ( 1 + 善) 婴：，7 ( f 一1 ) ( 3 7 ) ( 3 - 8 ) 3 1 2 二阶插值边界元法离散积分公式 边界积分方程的数值离散化格式有多种，奉文采用的是加权余量法对边界积分 方程进行离散化的伽辽金格式。将求解域的边界分割成有限个单元，然后在单元内 采用上述插值方法对电位函数妒和电位法向导数五= a c a a n 进行插值近似为 妒= m 纺1 五：圭m _ - 9 其中e 表示二阶插值单元，表示单元中结点的局部编号，n ，是插值单元的形状函 数。 由于采用伽辽金格式，所以取权函数为m ，将插值函数和权函数代入积分方程 式( 2 - 1 3 ) 的加权余量方程，可得 2 昭i f m 竹豳= i f l 鲁垃d s + “ “ “ ( 3 - 1 0 ) 莓i 莩i ，s 莩m 乃b r d s + ；i i 莩i s ；_ 纺! 警d s d s 式中，表示场点所在的边界面单元；d 表示源点所在的体单元( 等式右边第一项) 或边界面单元( 等式右边第二，三项) ；r 为源点到场点的距离：c o s o 表示场点到 源点的向量与源点处边界外法线方向夹角的余弦。 这是一组关于边界面上结点电位和电位法向导数的线性代数方程组。在一定的 边界条件下可以求解此方程组。 华北电力大学硕十学位论文 3 1 3 二阶插值边界元算法 由上述可知，选择某二阶插值单元e 就会有相廊的单元形状函数j 。设单元内 各结点的整体坐标为( 薯，片，弓) ，其中i 为单元内结点的局部编号。由于单元内任一点 的整体坐标可表示为局部坐标的函数，即 j = m 五，y = f 片。：= j ( 3 一“) fff 整体坐标对局部坐标求偏导，可得 言= 。a 。s v ；, 薯，嚣= 莓器乃，妻= 莩等毛 高2 7 z a 却, v , ，考；莩筹片，高= 莩筹刁 ( 3 一1 2 ) 二豢黼鞴锲矧_ 睁m= 南( ( 嚣荔一毒嚣 巳+ ( 责骞一毒骞) 勺+ ( 妻考一詈嵩) 巳 整体坐标面积元与局部坐标面积元的关系为 d s ：i , s l d 一山 华北电力人学硕+ 学位论文 5 2 3 实例分析三 给主体半径为0 2 m ，截面半径为0 0 4 m 的导体圆环加l o o v 电压，求圆环表面 的电场强度。奉文采用二维轴对称有限元方法的计算结果作为判断边界元方法计算 精度的依据，其中有限元单元数为8 1 7 4 ，结点数为1 6 6 9 0 ，计算所得的圆环血电场 强度幅值分布图如上图5 7 所示。 5 2 3 1 二阶插值边界元算法 ( a ) 单无数为1 9 2 ，结点数为3 8 4 的三角彤二阶插值单元 ( b ) 单元数为4 8 ，结点数为1 4 4 的l ! ! i 边彤二阶插值单元 图5 - 8 二阶插值边界元算法结果 一3 8 华北电力大学硕十学位论文 5 2 3 2 圆环面曲边四边形边界元算法 由于所求的是导体圆环模型，所以采用的是基于圆环坐标的曲边四边形边界元 法。其中，单元数为1 4 4 ，结点数为1 4 4 。用已经编好的计算机程序求得圆环面上 各点的电场强度如图5 - 9 所示。 图5 - 9 圆环面曲边四边形边界元算法结果 5 2 3 3结果比较 从图5 - 8 可以看出，( b ) 图的电场分布比( a ) 图的要匀称。在有限元法结果 中，电场强度值分布为( 3 3 8 6 5 ，9 3 5 3 7 6 ) ；在三角形二阶插值边界元算法结果 中，电场强度值分布为( 2 8 6 5 0 7 ，1 0 0 6 ) ；在四边形二阶插值边界元算法结果中， 电场强度值分布为( 3 0 7 2 9 8 ，9 0 8 9 8 ) ，显然，( b ) 图的电场强度整体分布更趋 近于高精度有限元计算的圆环面电场强度分布。下面只比较四边形二阶插值边界元 算法结果和圆环血曲边四边形边界元算法结果。如图5 - 1 0 所示。 从图5 - 1 0 看到，四边形二阶插值边界元和圆环面曲边四边形边界元这两种算 法的结果随着极角的总体分布趋势与有限元法是一致的，而且，两种算法的结果曲 线基本与有限元结果曲线能重合。相对而言，圆环面曲边四边形边界元结果曲线平 滑性较好些。上下极点处，四边形二阶插值边界元算法的误差分别为2 8 2 和 9 2 5 8 ：圆环面曲边四边形边界元算法的误差分别为i 4 和6 1 2 2 。 华北电力人学硕十学化论文 5 2 4 实例分析四 极角 图5 一1 0 圆环面电场强度随极角的分御 图5 一1 1 高精度有限元计算的圆环面电场强度 椎 华北电力人学硕十学位论文 主体半径为0 2 m ，截面半径为0 0 4 m ，相距为0 1 5 m 的两个导体圆环，上圆环 电位l o o v ，下圆环电位一l o o v 。用边界元法计算圆环面电场强度。 根据两导体圆环的形状和相对位置可知，上下导体表面电场强度应具有反对称 性，即上导体表面电场强度垂直于球面方向朝外，下导体表面电场强度垂直于球面 朝里，上下圆环对称点上电场强度大小相等。所以，本文用镜像原理，只对一个圆 环进行建模，再用相应的计算机程序算得圆环畸上备点的电场强度。为了评价边界 元计算结果的精度，根据轴对称性可以用二维轴对称电场计算导体圆环面的电场强 度。图5 一1 1 是用二维轴对称高精度有限元法( 考虑无限远处电位为零，削分单元 数为3 9 4 8 2 ，结点数为7 9 6 2 5 ) 计算所得的圆环面电场强度幅值分布图。 5 2 4 1二阶插值边界元算法 图5 - 1 2 二阶插值边界元算法绵粜 由于四边形二阶插值单元能更好地拟合边界曲面，从前面的实例分析可知三角 形单元的计算精度较低，所以在奉实例分析中只列出了单元数为4 8 ，结点数为1 4 4 华北电力人学硕士学位论文 的四边形二阶插值边界元的计算结果，如图5 一1 2 所示。 5 2 4 2 圆环面曲边四边形边界元算法 由于所求的是导体圆环模型，所以采用的是基于圆环坐标的曲边四边形边界元 法。其中，单元数为1 4 4 ，结点数为1 4 4 。用已经编好的计算机程序求得圆环面上 各点的电场强度如图5 1 3 所示。 图5 - 1 3 圆环面曲边四边形边界元算法结果 5 2 4 3结果比较 从图5 1 2 和图5 一1 3 可看到。圆环面曲边四边形边界元i ；i 。算得到的电场强度分 布很匀称，对称性很强，而四边形二阶插值边界元计算的电场强度在圆环底部分布 不是匀称的。 华北电力人学硕十学化论文 0 有限元结果 _ 圆环曲边四边形边界元结果 元结果 - 5 0 0 5 01 0 01 5 02 0 02 5 03 0 03 5 04 0 0 极角 图5 一1 4 圆环面电场强度随极角的分布 从图5 一1 4 看到，四边形二阶插值边界元和圆环面曲边四边形边界元这两种算 法的结果随着极角的总体分布趋势与有限元法是一致的，而且，两种算法的结果曲 线基本与有限元结果曲线能晕合。相对而言，圆环面曲边四边形边界元结果曲线平 滑性较好些。 5 3 本章小结 奉章首先介绍了三维静电场边界元程序的实现过程，二阶插值单元是奉文晕点 研究的对象。随后通过四个实例分析，针对不同的边界面形状，采取不同边界元算 法，比较了二阶插值边界元算法和球面曲边四边形边界元算法以及圆环面曲边四边 形边界元算法的计算精度。我们可以得出四边形二阶插值边界元是精度较高的一种 算法。同时，球面和圆环血曲边四边形边界元算法的计算精度在某种程度上要高于 四边形二阶插值边界元。 咖 啪 伽 瑚 啪 啪 晷伽 一e ；j | 华北电力人学硕七学位论文 第六章结论 本论文在总结电场计算方法的基础上，选择用边界元法求解三维静电场的开域 问题。由于三维静电场边界元法通常采用近似的数值积分法，要恕获得较为精确的 数值解就需要我们去研究和总结提高汁算精度的算法。 为了方便，奉文采用的是等参数单元。插值单元阶次的选取决定了有限单 元上逼近的精度，从而决定了整个问题的结果精度。选择合适的单元和单元内合 适的插值函数是提高边界元法计算精度的一个晕要课题。 本文深入研究三维静。场边界元法二阶单元的积分精度问题。分析了影响该积 分精度的相关因素，即单元形状、大小、个数和高斯点个数等对积分精度的影响。 并进而提出了解决奇异积分和改善积分精度的方法。编制了三维静电场边界元法的 二阶等参元有限元单元的计算软件。 本文还介绍了基于球坐标三维静电场曲边四边形边界元方法和基于圆环坐标 的三维静电场曲边四边形边界元方法。 通过对不同电极表面和导体表面采用线性单元、二阶单元和基于球坐标的曲边 四边形单元及基于圆环坐标的曲边四边形单元进行不同的嘲格剖分，奉文对比分析 了这四种边界元算法对计算精度的影响。 基于理论和计算结果，本文主要的结论如下： 1 ) 在等参数单元中，数值逼近的精度是用单元形函数来表达的，单元形函数 实际上包含了函数逼近和形状逼近两个方面，主要决定了面积单元上逼近的精 度。单元形函数的精度决定于插值函数和小单元的几何特性。例如，平面或曲 面边界等。单元上的插值函数的阶次取决于单元结点的个数。 2 ) 肖有曲面边界和复杂形状边界时，在函数逼近和形状逼近方面，由于采用 线性插值不够，所以针对研究的区域就要采用更高阶插值函数来逼近。本文考 虑到计算机内存大小，计算机速度等外在冈素而选择二阶单元，从实例结果可 得，二阶单元的计算精度的确要高于线性单元；四边形二阶单元的讣算精度要 高于三角形二阶单元。 3 ) 基于球坐标的曲边四边形边界元和基于圆环坐标的曲边四边形边界元这两 种算法中由于采用的是实际曲面的曲边单元和插值函数来满足函数逼近和形 状逼近这两方面，所以，它们的汁算精度要高于二阶单元。 4 ) 从实例分析中，我们除了能看到上述几种算法的计算精度差别之外，在计 算速度上，二阶单元由于增加了结点数量，所以其计算速度并不是很占优势。 华北电力人学硕十学化论文 1 何红雨 2 0 0 4 2 张丽萍 1 9 9 8 参考文献 电磁场数值计算法与m a t l a b 实现武汉：华中科技大学出版社 袁建生电磁场数值计算国际会议综述电t 技术杂志。第1 1 期 3 谢德馨，姚缨英，白保东，李锦彪三维涡流场的有限元分析北京：机械 工业出版社。2 0 0 1 4 3r w c l o u g h t h ef i n i t ee l e m e n tm e t h o do fs t r u c t u r a la n a l y s i s p r o c 2 n d c o n f e l e c t r o n i cc o m p u t a t i o n a s c e 。p i t t s b u r g h ，p a ，s e p t 1 9 6 0 5 w i n s l o w ，a m m a g n e t i cf i e l dc a l c u l a t i o ni na ni r r e g u l a rt r i a n g u l a rm e s h l a w r e n c er a d i a t i o nl a b o r a t o r y ( l i v e r m o r e ，c a l i f o r n i a ) ，u c r l 一7 7 8 4 一t ， r e v 1 1 9 6 5 6 p p 席尔维斯特，r