（电工理论与新技术专业论文）无单元方法在三维电磁场中的应用研究.pdf

屯 黔 1 乡 1 慧 声明尸刚 本人郑重声明：此处所提交的硕士学位论文无单元方法在三维电磁场中的 应用研究，是本人在华北电力大学攻读硕士于, - - - - z - , 位期间，在导师指导下进行的研 究工作和取得的研究成果。据本人所知，除了文中特别加以标注和致谢之处外， 论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得华北电力大 学或其他教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究 所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 学位论文作者签名：窒圭氐日期： 型：仝 关于学位论文使用授权的说明 本人完全了解华北电力大学有关保留、使用学位论文的规定，即：学校有 权保管、并向有关部门送交学位论文的原件与复印件；学校可以采用影印、缩 印或其它复制手段复制并保存学位论文；学校可允许学位论文被查阅或借阅； 学校可以学术交流为目的，复制赠送和交换学位论文：同意学校可以用不同 方式在不同媒体上发表、传播学位论文的全部或部分内容。 ( 涉密的学位论文在解密后遵守此规定) 作者签名：导师签名： 日期：蝼2 早益二。f b 华：i 匕i l 三力火! 学硕j = 学位论文 第一章绪论 无单元法是近几年发展起来的与有限元法相似的一种数值计算方法，形式上相 对简单、明确。在方法分类上，无单元法属于基于概率模型的计算方法，也是一 种偏微分方程弱形式的近似解法。它采用滑动最小二乘法构造形函数，从能量泛函 的弱变分形式中得到控制方程，并用拉格朗日乘子法或罚函数法等方法满足本征条 件，从而得到偏微分方程的数值解。 1 1 无单元法的研究概况 无单元法最早起源于力学领域。国际上对无单元法的研究兴起于上世纪9 0 年 代初，最早专门论述无单元法的文献【2 3 发表于1 9 9 4 年，从那时至今，研究热潮持 续不断。两份权威国际学术期刊c o m p u t e rm e t h o d s i na p p l i e dm e c h a n i c sa n d e n g i n e e r i n g 4 】和c o m p u t a t i o n a lm e c h a n i c s f 1 还为此出了专刊，可见国际学术界对包 括无单元法在内的无网格方法的重视程度。美国谣北大学( n o r t h w e s t e r nu n i v e r s i t y ) 的t b e l y t s c h k o ，w k l i u ，y k r o n g a u z ，d o r g a n ，p k r y s l 等几位学者在理论和 应用上对无单元法的发展做出了许多创新性的贡献。 无单元法的理论发展历经滑动最小二乘法【6 】( m l s ，m o v i n gl e a s ts q u a r e m e t h o d ) 到发散单元法【7 】( d e m ，d i f f u s ee l e m e n tm e t h o d ) ，再到无单元伽辽金法 ( e f g m ，e l e m e n t f r e eg a l e r k i nm e t h o d ，即无单元法) ，再发展到单位分解法( p u ， p a r t i t i o no f u n i t y ) 的主线，而最早的无网格法可追溯至用于研究天文现象的光滑粒 子法【8j ( s p h ，s m o o t h i n gp a r t i c l eh y d r o d y n a m i c s ) 。 无单元法的研究历史可大致概括如下： 它的前身是l u c y 在1 9 7 7 年提出的“s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ”方法( 光 滑质点流体力学方法，简称s p h 法) ，该法用于天体物理领域，并取得成功。 随后，在2 0 世纪8 0 年代，m o n a g h a n 等人发展了s p h 法，将其解释为核函数 法，并用来模拟流场中的激波强间断现象。 2 0 世纪9 0 年代，s w e g l e ，d y k a 等人提出了s p h 法不稳定的起因及稳定化方案： j o n h s o n 和b e i s s e l 等人也提出了一些对s p h 法的应变计算进行改善的方法。我国学 者对s p h 法进行研究的主要成果有国防科技大学贝新源、岳宗五等将s p h 法用于 高速碰撞问题等。 1 9 9 2 年n a y r o l e s 等人将滑动最小二乘法( m o v i n gl e a s t - s q u a r e sm e t h o d ，简称 m l s m ) 原理引入，提出了一种新的方法“d i f f u s ea p p r o x i m a t i o na n dd i f f u s e 华北电力人! 学预? h 学位论文 e l e m e n tm e t h o d ”( 散射元法，简称d e m ) ，山此无单元方法得到了进一步发展。 此后，b e l y t s c h k o 等人在1 9 9 4 年对d e m 进行了改进，基于滑动最小二乘法近似 和伽辽金( g a l e r k i n ) 原理，利用背景网格进行积分计算，提出了一种新的数值方 法“t h ee l e m e n t f r e eg a l e r k i nm e t h o d ”( 无单元g a l e r k i n 法，简称e f g 法) ，这种 方法采用滑动最小二乘法进行场的近似表达，虽然比s p h 计算费用高，但是具有较 好的协调性及稳定性。 1 9 9 5 年，美籍华人计算中心的学者w k l i u 教授根据积分再生函数的思想，提 出了“r e p r o d u c i n gk e r n e lp a r t i c l em e t h o d ”( 再生核质点方法，简称r k p m ) ，同时结 合小波分析中的伸缩尺度平移、多分辨率等特点，构造了“m u l t is c a l er e p r o d u c i n g k e r n e lp a r t i c l em e t h o d ”( 多尺度再生核粒子法，简称m r k p m ) 和“w a v e l e tp a r t i c l e m e t h o d ”( 小波质点法，简称w p m ) ，并实现了该方法的自适应分析。 同年，美国t e x a s 大学的著名学者o d e n , 年l l 他的学生d u a r t e 提出了h p c l o u d s 无网 格方法。后来，o d e n ，d u a r t e ，z i e n k i e w i c z 等又提出了“n e wc l o u d s b a s e dh pf e m ”。 这种方法借助于有限元网格，将其形函数作为单位分解函数，虽然该方法破坏了“无 网格”的特性，但由此带来了一些新的特性，为解决问题提供了方便。与h p c l o u d s 无网格方法相近的“h p m e s h l e s sc l o u d sm e t h o d ”是由波兰学者l i s z k a 等人提出的。它 最主要的特点是采用配点格式，不需要g a l e r k i n 格式中用于积分计算的背景网格， 是一种完全的无单元方法。 随着小波应用范围的扩展，2 0 世纪9 0 年代许多学者将小波用于场量的近似， 从而提出了“w a v e l e t g a l e r k i nm e t h o d ”( 小波伽辽金法，简称w g m ) 。 最近，著名力学学者a t l u r i 等人也对无网格方法进行了大量研究，提出了“l o c a l b o u n d a r yi n t e g r a le q u a t i o nm e t h o d ”( 局部边界积分法，简称l b i e 法) s o “m e s h l e s s l o c a lp e t r o v g a l e r k i nm e t h o d ”( 简称m l p g 法) 【9 1 。这两种方法都是基于滑动最小二 乘法原理对场函数进行近似，并且在积分时不需要背景网格，是完全的无单元计算 方法。 国内对无单元法研究是从1 9 9 8 年清华大学的周维恒、寇晓东用无单元法解决 结构开裂问题【1 0 ，1 1 1 开始的，2 0 0 4 年，胡云进、周维恒、寇晓东借鉴二维无单元法的 理论和实现技术，对基于正交基函数的三维滑动最小二乘法、已知位移边界条件的 处理以及三维无单元法基本方程的推导和实现方法进行了讨论。上述成果都是应用 于力学领域。 无单元法在电磁场领域的发展相对滞后。 1 9 9 3 年，法国学者y m a r e c h a l ，j l c o u l o m b ，g m e u n i e r ，g t o u z o t 应用扩 散单元法( d i f f u s ee l e m e n tm e t h o d ) 求解电磁场问题。这可以看作是无单元思想在 华北l u 力人学顶j 卜：学位论文 电磁场领域的首次尝试。 1 9 9 8 年，日本学者v l a t k oc i n g o s k i 等人将无单元伽辽金法( e f g m ) 首次应 用于电磁场计算领域，用于解决电压互感器中的二维静电场问题。2 0 0 0 年，上述的 日本学者将e f g m 和有限元法( f e m ) 结合，分析了二维电磁场问题。但用e f g m 与 f e m 耦合法处理强加边界条件和介质交界面条件时，虽然保证了系数矩阵的对称正 定性，但在交界区的形函数过于复杂。 1 9 9 9 年，巴西的v i a n a 等人用拉格朗日( l a g r a n g e ) 乘子法对强加边界条件进行 了处理，并给出了二维磁场的算例。但l a g r a n g e 乘子法最大的缺点是它引入了新的 未知量，使得离散方程的系数矩阵不再是正定的，也不是带状的，增加了方程求解 的时间和难度。 1 9 9 9 年，法国的c h e r a u l t 等人对无单元伽辽金法( e f g m ) 的边界条件和介质 交界面条件处理问题进行了初步研究，给出了处理强加边界条件和介质交界面条件 的几种处理方法，如l a g r a n g e 乘子法、跳跃函数法等，但文中并没有给出具体算例， 而且给出的仅仅是一维跳跃函数。 1 9 9 9 年，浙江大学杨仕友、倪光正在分析国内外研究现状，综合前人工作的基 础上，将一般紧支集正交小波应用于电磁场的数值分析和计算中，提出了小波一伽 辽金法中任意点关联系数值精确的时域计算方法。 1 9 9 7 年，河北工业大学刘索贞，杨庆新等人借鉴无单元g a l e r k i n 法在力学场应 用的成功经验，开始对其在电磁：场领域的应用进行了研究，刘素贞在其硕士论文中 研究了二维静态电场磁场的无单元方法，进行了简单的工程实例计算，如平行板电 容器中的电场、长直接地金属槽中的电场、电机气隙和槽部磁场、变压器中二维磁 场、长直矩形载流导线的磁场。研究工作取得了初步的进展，但是距离工程实际应 用还有很大差距。 此外，最近还有学者用无单元g a l e r k i n 法解决冲击涡电流问题。国内外学者在 应用核函数的无单元法解决电磁场问题上也取得了一些成果。 迄今为止，无单元方法在电磁场领域取得一些研究成果，但该方法在电磁场领 域的研究仍处于起步阶段，还存在着许多亟待解决的问题。 1 2 本课题的研究意义和目的 有限元法等传统方法是求解电磁场边值问题强有力的数值计算方法。这类方法 多数以单元作为基本概念，将场域空间离散成单元的形式，采用剖分插值的方法近 似表征场域函数。对于单元的形状，有一定的要求，不能出现过大的钝角和过小的 锐角。在场函数变化剧烈或者介质参数变化剧烈的位黄，必须进行相对比较密集的 毕：i 匕电力人学硕j ：学位论文 剖分，费时费力。尤其在处理三维问题时，网格剖分成为前处理阶段最为繁琐的工 作。当这类方法用于自适应计算或动态模型分析时，一般要不断更新网格。虽然目 前已出现一些网格生成器，但在准备数据阶段还是占用机时多，不方便。 与传统的有限元法相比，无单元法具有以下一些显著的优点： ( 1 ) 只需节点信息，不需将节点连成单元。要求数据简单，完全不需要事先 对计算区域进行网格剖分，在计算过程中可方便地增删节点，极大地简化了前处理 工作。 ( 2 ) 由于采用滑动最小二乘法构造形函数，故所得的场函数光滑性好，计算 精度高，可以满足连续可导甚至高阶连续可导。对于要求离散精度较高的位置可以 方便地加密节点，而不用顾虑单元形状的问题。 ( 3 ) 在理论和应用的各个环节上均具有较高的开放性，便于进行控制和扩展。 无单元法中基函数、权函数、影响半径等的选取具有较大的随意性，对边界条件的 处理方法也有多种。基于这种开放性，无单元法的前景是非常广阔的，特别是一些 前沿科学的引入( 如小波分析) ，必将使无单元法在理论和应用上取得质的飞跃。 正是由于具有这些显著优点及巨大的发展潜能，无单元法得到了日益众多的关 注，并取得了一定的发展。在电磁场领域的研究成果很多，但主要集中在解决二维 电磁场问题，在三维电磁场中的应用很少。 1 3 本文的主要研究内容 本文的主要研究内容有以下几个方面： ( 1 ) 借鉴无单元法在二维电磁场中的基本公式和实现技术，从伽辽金弱变分 原理出发，推导无单元法在三维静态电磁场中的计算格式，分别给出拉氏乘子法和 罚函数法处理边界条件的一般公式； ( 2 ) 探讨无单元法边界条件及介质交界面的处理方法； ( 3 ) 通过分析算例，研究影响无单元法计算精度的主要因素，包括基函数、 权函数、影响半径和节点分布密度： ( 4 ) 在m a t l a b 软件环境下，实现三维静态电磁场的无单元法的计算程序， 并将结果与有限元的计算结果进行比较。 华：i 匕电力人：学硕十学位论文 第二章无单元法及其在电磁场计算中的应用 本章借助无单元法在二维静电场中的基本公式和实现技术，将之推广到三维静 电场数值计算中。首先介绍理论基础一滑动最小二乘法的基本原理，再从伽辽金 弱变分原理出发，推导出无单元法在三维电场和磁场下的计算格式，另外也涉及了 权函数、影响半径的选取方法和边界条件的处理方法等内容。 2 1 理论基础 2 1 1 滑动最小二乘法 滑动最小二乘法( m o v i n g l e a s ts q u a r e ，简记：m l s ) 是无单元法的数学基础， 其基本原理是通过若干个节点上的值，拟合出一个函数。不同于一般最小二乘法采 用全局最小二乘逼近，滑动最小二乘法借助一个影n 向域有限、单调递减的权函数， 采取了局部最小二乘法逼近的策略，在各个局部均强调了数值解与真实解之间的拟 合关系，它的待定系数是空间坐标的函数，因此可以构造出更加精确的函数形式。 在空间场域q 中，设场函数为u ( x ) ，对于二维空间，x = ( x ，) ，) t ；对于三维空间， x = ( x ，y ，z ) t 。若给定u ( x ) 在f 个节点上的值为 “+ ( x i ) = 1 1 7 ， i = 1 ，1 ，咒 ( 2 一l 一1 ) 则以这n 个节点上的值，应用滑动最小二乘法可拟合出场函数的个近似函数( 也 称试函数) u 6 ( x ) ： ( 2 1 2 ) 式中p t ( x ) 是m 维基函数，一般选择多项式，当然，也可以使用其他函数作为基函 数。例如在分析奇异性问题时，可以将奇异函数作为一个基函数。a ( x ) 是m 维线性 系数矩阵，它是关于场点( 工，) ，) 或( r ，j ，z ) 的函数 a ( x ) t = k 。( x ) a 2 ( x ) a 3 ( x ) 】 对于一维情况： 线性基p t ( x ) = 1 ，工 ，7 2 = 2 ( 2 一l 一3 ) ( 2 一l 一4 a ) x ，i i 、 a 、， x ，_、 t p l i 、l ， x i 以 、， x ，i p 川削 一一 、， x ，l “ 华北电力人! 学顸7 1 7 ”j - 位论文 二次基p t ( x ) = 1 ，x ，x 2 ，m = 3 三次基p t ( x ) = 1 ，x ，x 2 ，x3 ，m = 4 对于二维情况： 线性基 p t ( x ) = 1 ，x ，y ，2 = 3 二次基 p t ( x ) = 1 ，x ，y ，x 2x y ，y 2 】，m = 6 ( 2 1 4 b ) ( 2 1 4 6 ) ( 2 1 4 e ) 三次基 p t ( x ) = 1 ，x ，y ，x 2x y ，y 2 , x 3 ，x 2 y ，砂2 ，y 3 ，7 2 _ 1 0 ( 2 1 4 f ) 对于三维情况： 线性基 p t ( x ) = 1 ，x ，y ，z ，m = 4 ( 2 一l 一4 9 ) 二次基 p t ( x ) = 1 ，x ，y ，z ，x 2 , y 2 , z 2x y ，x z ，y z ，m = 1 0 ( 2 1 4 h ) 三次基 p t ( x ) = 1 ，x ，y ，z ，x 2y 2 , z 2 , x y ，x z ，y z ，x 3 ，y 3 , z 3 , x 2 y ，x 2 z ，y 2 x ，y 2 z ，z 2 x ，z 2 y ，x y z ， m = 2 0( 2 1 4 i ) 式( 2 。1 2 ) 是全局近似，而在近似计算中一般不采用全局近似，一方面收敛性不 容易控制，另一方面，也不利减小计算量，所以要进行分片拟合，即要进行局部近 似。其对应的局部近似表达式为： “6 ( x ，x ，) = p 小加小) = p t ( x ，) a ( x ) ( 2 1 5 ) 为了确定a ( x ) ，在局部范围内构造带权重的上2 范数j ( x ) ，并使j ( x ) 取最小值： j ( x ) = w ( x - x ，) “( x ，x ，) - u + ( x ，) 2 ， = w ( x x ) p t ( x ，) a ( x ) - u + ( x ，) 2 ( 2 1 6 ) 式中x ，是x 影响域( s u p p o r to f i n f l u e n c ed o m a i n ) p 勺的点，共有咒个。这门个节点不必 连成单元，这是无单元法有别于有限元法之处；w ( x b x l ) 是权函数，其作用是控制 距离场点不同的节点对场点影响的大小，对影响域外的x ，w ( x x ，) = 0 ，对影响 域内的点x ，w ( x - x t ) 0 ；t l * ( x ，) 是“+ 在节点x ，处的值。 将式( 2 1 6 ) 写成矩阵形式 j ( x ) = ( p a u 。) 一w ( x ) ( p a u4 ) ( 2 一l 一7 ) 华：i 匕【匕力人学顺十：学位论文 式中 p = w ( x ) = p l ( x 1 ) p a ( x 1 ) p 川( x 1 ) p 1 ( x 2 ) p 2 ( x 2 ) p 川( x 2 ) p l ( x ，：) p 2 ( x ，) p 川( x ，) w ( x x 1 ) 0 0 0 w ( x x 2 ) 0 00 w ( x x ，1 ) d * ( x 1 ) z f + ( x 2 ) 当j ) 取最小值时，由梨：o ，得到 0 a i x ) a ( x ) a ( x ) 一b ( x ) u + = 0 式中，中间变量a ( x ) 和b ( x ) 分别为 a ( x ) = p x w ( x ) p b ( x ) = p t w ( x ) 当a 。1 ( x ) 存在时，解出 a ( x ) = a 。1 ( x ) b ( x ) u + ( 2 一l 一7 a ) ( 2 1 7 c ) ( 2 1 8 ) ( 2 1 8 a ) ( 2 1 8 b ) ( 2 1 9 ) 每个局部近似式( 2 1 5 ) 的a ( x ) 都对全局近似式( 2 1 2 ) 有贡献，将式( 2 1 8 a ) 矛n ( 2 1 ：8 b ) 代入( 2 1 2 ) 中，有 “( x ) = p t ( x ) a 一( x ) b ( x ) u + ( 2 1 1 0 ) 注意一点，m l s 不满足克罗内克6 条t ：- 1 ：( k r o n e c k e rd e l t ac r i t e r i o n ) ，即不满足 插值条件。由于采用滑动最小二乘法拟合，所以在x ，处，由式( 2 一l 一1 0 ) 得到的 “( x ，) u ( x ，) 。因此，在求解电磁场边值i 口- j 题时，无单元法对第一类边界条件的处 理将有别于有限元法，而是采用拉格朗r 乘子法或罚函数法等的特殊方法。 l p ：i ll u 力人学硕：卜学位论文 2 1 2 形函数及其物理意义 可将式( 2 1 1 0 ) 写成下面形式 “6 ( x ) = n ( x ) u + 式中n ( x ) 称为形函数，其定义为 ( 2 1 1 1 ) n ( x ) = p t ( x ) a 一( x ) b ( x ) = 。：，】 ( 2 1 1 2 ) 尽管式( 2 1 1 1 ) 是定义在所有节点上的和，但由于b ( x ) 是一个稀疏矩阵，大多数的m 的值都是零；只有在影n 向域内节点的形函数不为零。 形函数n ( x ) 的微分形式为 n ( x ) ，= ip t ( x ) a 。1 ( x ) b ( x ) ii = p t ( x ) ，a _ ( x ) b ( x ) + p t ( x ) a - 1 ( x ) ，b ( x 、fp t ( x ) a 一( x ) 【b ( x ) 】， l 广1 一 r1 ( 2 一l 一13 ) 为了计算式中a 一( x ) f ，由单位矩阵可以导出a 叫( x ) ，的表达式 i = a a j i ，f = ( a a 。) ， j 0 = a f a 一+ a a i ： ( 2 1 1 4 ) j a 叫( x 1 f = 一a 一( x ) a ( x ) f a 。1 ( x ) 式中f 代表空间坐标变量x 中的元素x 、y 或z 。 形函数的物理意义：形函数n ( x ) 可分为三个独立的部分，即p t ( x ) ，a - 1 ( x ) 和 b ( x ) 。其中p t ( x ) 表示插值点的位置信息；a 。1 ( x ) 是以插值点为圆心的区域内( 即影 响域) ，所有插值基点的坐标各自进行包括加权在内的运算后叠加而成的；b ( x ) 中 的每一列可看成是表示对应的插值基点与插值点x 之间的分配关系。因此，无单元 法的形函数可以理解为，即对于插值点x ，先由其影响域内的所有插值基点所形成 代表整个影响域的a ( x ) ，并求逆，再利用b ( x ) 将这一整体信息分配到域内的各插值 基点，就可得到插值点与这些基点问的插值函数值。 2 1 3 权函数 权函数是无单元法的精华所在，它对场函数的光滑性，计算量的大小，计算结 果收敛的快慢都有较大影响。若选取得当，可使近似解和精确解在某一划分的区间 内更好的拟合，所以权函数的选取非常重要。目前，权函数的选取还没有一种通用 的方法，主要依靠使用者实践经验的积累，但其选取仍需满足以下基本原则1 ： 华：i 匕l 也力人学硕十学位论文 ( 1 ) 权函数应非负； ( 2 ) 权函数应连续可导，以保证近似场函数连续可导。 ( 3 ) 权函数应单调递减，某节点的； 叉函数在该节点取得最大值，由近及远逐 渐衰减，且在影响半径之外值为0 。 权函数的非负性，保证了待求场点周围节点对该点的贡献为正；连续可导性能 够影响场函数的特性，如权函数有k 阶连续偏导，则形函数也将存在k 阶偏导，也 就是说，只要权函数选择的适当，场函数可以获得任意阶连续偏导；单调递减性是 指在影响半径之外权函数为0 ，这样保证了在求解a ( x ) 时a 一( x ) 的存在性。 2 1 3 1 常用权函数形式及比较 满足以上基本规则的前提下，无单元法常用的权函数类型有以下几种【2 j 1 1 2 、1 3 】： 指数型( e x p o n e n t i a l i ) ，、f e 一( 口) j 1 以d31 o刚 负指数型( e x p 。n e n t i a l i i ) “s ) = 百玎 f e 一( 嬲) 。一e 一口” 0 ， 三次样条( c u b i cs p l i n e ) 四次样条( q u a r t i cs p l i n e ) 五次样条( q u i n t i cs p l i n e ) ，埙s，=43兰345-+4s2芝s43j3，s 以班1 - 6 s + p s 4 吣， 1 - 1 0 s + p 6 5 5 有理式形式( r a t i 。n a lf 。r m u l a ) “s ) = ( 1 k g s 2 + 乒) 圆锥型权函数( c o n i c a l ) s 1 s 1 j 1 s 1 s 1 2 1 2 1 s 1 s 1 s 1 s 1 s l s 1 式中s = 0 r 1 ，其中0 = i lx x 川，l 。，为影响半径。式( 2 1 15 a ) 中口为控制系数； 式( 2 1 1 5 b ) q b 口为控制相对权重的参数；式( 2 1 1 5 f ) 中k 为正整数，孝为一很小的正 数；式( 2 一l 一1 5 9 ) 中k 为控制系数。 华北电力人! ! 浮硕十学位论文 常用权函数 图2 1 权函数图形 图2 1 给出了权函数的图形。通过分析可以看到，各种权函数的特性大体相同， 但是所具有的连续可导性不尽相同。三次、四次、五次样条函数型权函数和有理式 型权函数为严格意义的连续可导。而另外三种形式的函数，即指数型、负指数型权 函数和圆锥型权函数，在数学上都不是严格连续可导的，但是只要参数选取得合适， 可以保证近似连续可导，这对于数值计算结果的影响很小，可以忽略。对于一般精 度而言，参数选取范围可确定为：指数型，口不大于o 4 ，最佳为o 3 左右；负指数 型，2 5 t 6 ；有理式型k 可取一正整数，如2 ，4 ，6 等，f 可取小于或等于1 的 整数；圆锥型，固 1 式中；：坚型。其函数图形如图3 3 所示。 该跳跃函数能较好地模拟导数的 i 连续，但在二维区域罩面不具有紧支集特 华，i li l = ! 力人学硕? 卜学位论文 性，必须对其进行修正。在二维区域中，近似修- _ f 为 “( x ) = ，( x ) “，+ 6 ( s ) 甲( r ) = ，( 工) “，+ 钆m ( s ) 甲( r ) ( 3 1 3 ) llj 式中m ( j ) 是用弧长为自变量的一维形函数；轧疋，- 歇工l - t 一：。t l 跳跃强度的参数。 q 4 0 2 一夕 。i 一0 ；。一 6 j 0 1 5 j 。“i 。 0 2 ： - 0 4 ： 图3 3 一维跳跃函数 3 2 复杂边界的处理方法 图3 4 曲线边界近似 当求解区域出现比较复杂的边界时，会给网格剖分带来极大不便。为了不使方 程发生奇异，必须进行比较密集的剖分，从而大大增加了信息量，这是所有以单元 为基础的算法的弊端。无单元方法不需要节点构成单元，所以节点分布可以非常灵 活。在复杂的边界附近可以根据需要加密节点，只要保证不使相邻节点间距差别过 大即可保证计算精度。 另外，边界上没有必要布置过多的节点。原因是无单元法采用的是拟合方法而 不是有限元方法采用的插值方法，所以节点：倍置自由得多，只要保证边界的附近和 内部有足够多的节点分布即可，对节点的位置不做苛刻的要求。 关于第一类边界条件的处理，无单元方法是通过边界积分的方法加入的，对于 二维情况，做线积分，对于三维情况，做面积分。当遇到曲线边界时候，本文采用 折线代替曲线的方法( 三维为平面代替曲面) 。然后分别作各条线段( 平面) 的积 分，如图3 4 所示。 华一匕电力人学硕士学胁0 论文 3 3 算例 考虑一个简单的二维静磁场模型，其平面示意图如图3 5 所示。图中存在两利， 材料，一种是相对磁导率为1 5 的线性材料，“= 1 5 a t 。，另一种是相对磁导率为1 的空气，u ，= 肌，流过的面电流密度j = 1 5 0 0 a i l l 2 。以矢量磁位a 为求解量进行计 算，零磁势参考点本应该选无限远处，但为了简化计算，假定外沿的磁势为零，此 为第一类边界条件，该模型结构沿x 轴有齐次的第二类边界条件成立。对于这两种 介质交界面的问题，采用可视性准则的方法进行处理。 利用无单元法计算出各节点的a ，再由b = v a 求出b ，并绘制出该问题的磁 场分布，如图3 - 6 所示。 图3 5 二维磁场模型的计算区域 3 4 小结 图3 - 6 磁场分布图 本章主要介绍了处理介质交界面的三种方法，包括可视性准则、l a g r a n g e 乘子法 和跳跃函数法，阐述了各法的基本原理，分析了它们的优缺点。通过对二维静磁场模 型的计算，实现了用可视性准则的方法处理介质交界面问题，得到了比较满意的结果。 另外，对复杂边界的处理方法也做了简单的讨论。 华：l 匕电力人! ! 弘硕士学位论文 第四章对影晌计算精度的四个主要因素的研究 第二章中已对无单元法的基本原理作了系统的阐述。通过对计算原理的研究， 可以发现，影响无单元法计算精度的因素主要有以下几个方面： ( 1 ) 权函数 ( 2 ) 基函数 ( 3 ) 影响半径 ( 4 ) 节点分布密度 在这一章，将结合算例对以上四个影响无单元法计算精度的主要因素进行详细 讨论。 4 1 计算模型 如图4 1 所示，单位为米。根据分离变量法，槽内的电位的解析解为【1 8 】： “c x ，y ，= ：妻。s i n ；焉y 三竺兰耋擎 。4 一。， ：争s i n2 n s i n h2 n ：c ( 0 5 - ：r ) x 图4 1 长直导电槽截面 4 2 权函数对计算精度的影响 权函数在无单元法中起者至关重要| ，j 作用，是无蚺元法的核心，它的选取将直 27 白一 白 华：j 匕l 包力人：j o 硕一忙学位论文 接关系到无单元法的计算精度和计算的复杂性。在第二章中已对权函数的选取原则 及权函数的类型作了详细的阐述。这里分别选用第二章中式( 2 一l 一15 a ) 式( 2 一l l5 9 ) 作为权函数来计算图4 一l 中的算例，并将结果进行对比。 将场域用1 1 5 6 ( 3 4 x3 4 ) 个均匀节点离散，积分采用1 0 8 9 个矩形积分子域， 户四阶高斯积分，选择二次基作为基函数。采用上述七种权函数分别进行计算，权函 数参数选取 2 0 】分别为：有理式型( k = 4 ，f = 1 ) ：指数型( 口= o 3 5 ) ；负指数型 ( 口= 2 7 5 8 ) ；圆锥型( k = 1 ) 。将计算结果与解析解进行比较，表4 1 中列出计算 点的最大相对误差、均方差和相对c p u 时间。 表4 1 不同权函数下计算结果的比较 从表4 1 中的误差比较和c p u 时间比较可以看到，权函数的特性不同，计算的 精度差别较大。具有较好连续性和可导性的权函数，计算精度较高；反之，如圆锥 型权函数，就会造成较大误差；不同权函数形式下的计算时间略有差别，形式简单 的权函数的计算时间较短，如指数型权函数、圆锥型权函数，反之，形式复杂的函 数，计算速度比较慢，如三次样条函数。当计算的规模变大以后，这牙, f - i , - i - 算效率的 不同将会很明显。 总之，对于权函数的选取要进行比较全面的考虑，根据不同的场域特性确定不 同的权函数及其参数。本文通过大量的计算，总结出在电磁场分析中，有理式函数、 样条函数和指数型权函数均可以得到比较理想的结果。 4 3 基函数对计算精度的影响 在无单元法中，基函数的选用主要有线性基、二次基、三次基三种，基函数选 取的好坏对计算精度有显著的影n 阳。在同样的影响域、同样的节点密度和同样的权 函数情况下，选用不同的基函数分别进行计算。汁算结果如表4 - 2 所示。 2 s 华，i l l 也力人! 学硕? 卜! 学位沦文 表4 - 2 不同基函数卜汁算结果的比较 从表4 2 中可以看出j 选用三次基函数，计算精度比较高且计算效率偏低；选 用二次基函数，计算精度有所下降，但可提高计算效率；取线形基函数，计算结果 的误差要比前两个基函数大。一般情况下，取二次基函数就可以得到比较理想的结 果。 4 4 影响半径对计算精度的影响 无单元法是函数遐近，从这一角度看，影响半径的作用在于在影响域内包含足 够多又尽可能少的节点。“足够多”在于要有足够的采样点来满足场函数的连续性 要求并确定近似解中的待定系数；“尽可能少”在于要保证和突出邻近节点问的函 数关系免受其它节点的过多影响，从而使近似解更好地逼近真实解。因此需对影响 半径进行合理地选取。 选取具有一定代表意义的影n 向半径进行研究，对其计算结果进行对比，见表4 3 所示。 表4 3 不同影响半径的计算结果的比较 由表4 3 可以看出，在权函数和基函数选定及节点密度不变的情况下，当影响 半径较小时，计算结果与精确解相比误差较大，出现失真的情况。但这并不意味着 影n 向半径越大，计算结果越好。d ：l 表中可见，当影i 响域0 0 6 时，计算结果也会产 华北l l l 力人：学硕：卜：学位论文 生较大的误差，出现失真的情况。+ 从计算结果来看，当影n 向半径选为o 0 3 0 0 4 之 间时，计算结果较好。表中最后一行的影响半径是采用第二章提到的经验公式 ( 2 1 1 6 a ) 算出的，由此得出的汁算精度较高，计算效率也较高。 4 5 节点分布密度对计算精度的影响 节点分布密度对无单元法的计算精度有着比较明显的影响，为了研究节点分布 密度对计算精度的影响，基函数选用线性基，权函数选用指数型，影响半径选用经 验公式算出的值( 注：式( 2 1 16 a ) 中影向半径，j 7 ，与节点密度c 有关) 。节点数分 别选用1 2 1 ( 1 1 1 1 ) 、4 4 1 ( 2 1 2 1 ) 、9 6 1 ( 3l 3 1 ) 、16 8 1 ( 4 1 4 1 ) 四种情况进 行计算。如表4 4 所示，是四科，不同节点密度的计算结果的对比。 表4 4 不同：俸点密度的计算结果的比较 由表4 4 可以看出，在权函数及基函数选定的情况下，节点的分布密度对无单 元法的计算精度有着较大的影响，总体规律为随着节点密度的增大，计算精度有所 提高，但是当节点数过多，要以增加计算时间为代价，计算效率随之降低。如果在 影响半径不变的情况下，当节点数过少或过多，都有可能造成计算结果失真。 4 6 小结 本章对影响无单元法计算精度的四个主要因素权函数、基函数、影响半径和节 点分布密度进行了探讨，通过分析算例，得出了一些规律： ( 1 ) 在基函数选定，节点密度和影向半径不变的情况下，七种权函数都可以 满足计算精度的要求，各有优缺点； ( 2 ) 在权函数选定，节点密度和节点影响域不变的情况下，高次基函数比低 次基函数计算精度要高，但对计算结果影响不大； ( 3 ) 在权函数和基函数选定，节点密度不变的情况下：影响半径的作用在于 在影响域内包含足够多又尽可能少的：肖点，其值通常由经验公式求出； ( 4 ) 在权函数和基函数选定影响：仁径不变的情况下，节点分布密度对训：算 ： n 华北l 岜力人! 学颂i ：学位论文 结果影响较大，总体规律为随着节点密度的增大，汁算精度有所提高。 以上只是一些比较简单的规律，并不能单一的拿出来说明某个问题，这四个主 要影响因素是相互相关的、相互影响的。针对具体问题需要具体研究，不能为了追 求精度，而盲目的增加节点影n 向域和节点分布密度或者提高基函数的次数等等，这 有可能使计算结果失真。 华：f i = l 也力人! 学硕? = 学位论文 第五章计算实例 基于无单元法的基本原理，本章自行设计了三维电磁场无单元法计算程序，并应用 于分析三维静电场算例，所得结果与有限元法的计算结果进行比较，分析了计算结果的 误差，验证本算法的正确性和有效性。 5 1 计算模型 下面对一简单的三维静电场模型进行分析。图5 1 表示一方形电极和外屏蔽的 模型。根据对称性，只计算第一卦限内的部分，如图5 2 所示。图中单位为米。内 导体电极的电位取10 伏，外层屏蔽电位取0 伏，之间的区域认为是真空，无体电 荷分布。 图5 1 屏蔽电极模型 边界条件如下： 图5 2 计算区域及几何尺寸 平醯= 1 0 ，0 y 1 0 ，0 z 1 0 ：= o v 平面j ，= 1 0 ，0 x 1 0 , 0 z 1 0 ：矽= o v 平面z = 1 0 , 0 y 1 0 ，0 x 1 0 ：= o v 平面x = 5 , 0 y 5 , 0 z 5 ：= l o v 平面y = 5 , 0 x 5 , 0 z 5 ：= l o v 平面z = 5 , 0 y 5 , 0 x 5 ：矽= i o v 面x = 0 ，y 5 & z 5 ： 面y = 0 ，x 5 z 5 ： 却，、 2u a n a 矽，、 2u 功2 面z ：b ，再曩再石：掣：o o n ：3 2 华：i 匕电力人学硕十学位论文 将二维电磁场无单元方法进行推广，不难得到三维场的计算方法。基函数的选 取可以参照式( 2 1 4 9 ) ( 2 一l 一4 i ) ，其他公式进行相应的扩充即可。刺于三维静 电场问题，选取标量电位, 作为待求量。权函数的影响域变成了球形，对应的影响 节点个数明显增加。计算方程系数阵s 中的元素的积分变成了体积分，边界积分变 成了面积分。这样，积分点的个数也变得非常庞大，计算代价随之增大。 5 2 参数选取 ( 1 ) 基函数选取线性基，如式( 2 1 4 9 ) ；四阶高斯积分；3 0 3 2 个节点均匀分布。 ( 2 ) 权函数选取有理式形式( 2 1 1 5 f ) ，有理式函数为严格意义的连续可导函 数，计算精度较好。k 及s 的选择具有一定程度的任意性，但选取得好可提高计算 精度。有理式权函数具有k 一1 阶导函数连续，但是当k 取得较大时，影响域明显缩 小，容易造成计算失真，一般k 取值不超过4 js 为一个 圈、的正数，其值越小，权 函数的非零影响域越小。一般地，保证权函数在节点的取值应为1 ，s 可取为1 。对 于本文的模型，取k = 4 ，s = 1 0 时，计算结果误差最小，如表5 1 所示。 表5 1k 、s 在不同取值下的计算结果的误差 尼s 最大相对误差绝对值( )均方差( ) o 28 3 3 e 一0 43 8 7 e 0 4 o 57 9 5 e 一0 43 0 8 e 一0 4 2 o 87 6 8 e 一0 42 5 4 e 一0 4 1 07 4 0 e 一0 42 0 6 e 一0 4 o 28 1 1 e 一0 43 0 5 e 一0 4 4 o 57 8 4 e 一0 42 4 9 e 0 4 o 87 3 2 e 0 41 9 6 e 一0 4 1 o7 0 7 e 一0 41 8 2 e 一0 4 0 28 2 6 e 一0 43 1 7 e 一0 4 o 57 8 7 e 一0 42 8 4 e 一1 9 4 ， 6 o 87 4 2 e 0 42 0 1e 0 4 1 07 2 6 e 一0 41 8 3 e 一0 4 ( 3 ) 影响半径也需要适当选取，借鉴力学场的经验公式1 9 1 ，在节点均匀分 布时可取为 3 口，z 2 可丽 ( 5 1 ) 式中：n = 4 ( 线性基) ，10 ( 二次基) ，2 0 ( 三次基) ；c 为节点分稀密度；a 为大 于1 的系数，_ 般取4 6 ，本文取a = - 4 。- 1 i - + 点不均匀分布时，各点的，j ，需分别确 定，但原则同上式，只不过此时c 为变量。 华北l u 力火! 学硕+ 学位论文 ( 4 ) 在此例中采用罚函数法处理第一类边界条件。它的优点在于它保证了系 数矩阵的正定性和带状分布。理论上，罚因子越大，所得的结果越精确，但实际上 c 值取得越大计算起来越复杂，取得较小又不能满足精度要求，所以目前罚因子c 的选取尚无统一标准。这里取c = 1 0 0 0 。 5 3 程序结构 选取m a t l a b 作为无单元算法的开发环境，依据第二章的基本原理，可以将 无单元法的实现步骤总结如下： 生成节点并计算节点的影响半径。 确定高斯点位置及积分权。 对高斯点循环，扫描高斯点的影响节点，通过高斯积分，计算高斯点对整体平衡 方程的贡献，并将其集入整体平衡方程中。 如高斯点循环完毕，则执行，否则继续。 解整体平衡方程，得出节点电位值。 后处理，描绘等位线。 对应以上步骤，在m a t l a b 环境下，计算程序包括一个主程序( 脚本文件) m a i n m 和四个子程序( 函数分件) 包括：d a t ap r o m 准备数据一一布置节点，积分 点，计算积分权；g a u s s m 计算高斯积分； w e i g