（电机与电器专业论文）电磁场数值计算中的无网格方法研究.pdf

竺要翌三查主堡主兰堡笙兰 a b s t r a c t t h ef e mh a sb e e ne s t a b l i s h e da s a v e r yp o w e r f u l a n dm a t u r en u m e r i c a l t e c h n i q u e i np d e s c o m p u t a t i o n h o w e v e r ，t h e f e me l e m e n t g e n e r a t i o n c a n s o m e t i m e sb ev e r yl a b o r i o u sa n dt i m e c o n s u m i n g m o r e o v e r ，f o rp r o b l e m ss u c ha s a d a p t i v ec o m p u t a t i o n a n di n v e r s e s h a p eo p t i m i z a t i o np r o b l e m s ，f e mu s u a l l y r e q u i r e sr e m e s h i n g t oo v e r c o m et h e s ed i f f i c u l t i e s ，m e s h l e s sm e t h o d sh a v eb e e n d e v e l o p e db y r e s e a r c h e r si nm e c h a n i c a la n dm a t e r i a lf i e l d d u r i n g t h e p a s t t w o d e c a d e s i ti sb a s e do nn o d a ld i s t r i b u t i o na n dd o e sn o tn e e dm e s h g e n e r a t i o n e l e m e n t f r e eg a l e r k i nm e t h o d ( e f g ) i st h em o s tt y p i c a lm e s h l e s sm e t h o d s c o m p a r i n g w i t ho t h e rm e s h l e s s m e t h o d s ，i t s h o w sa t t r a c t i v e a d v a n t a g e s ( e g n u m e r i c a ls t a b i l i t y ，c o n v e n i e n tp o s t - t r e a t m e n t ，f a s tc o n v e r g e n c e ) i nt h i sp a p e r ，t h e b a s i cp r i n c i p l eo fe f ga n dt h em o v i n gl e a s ts q u a r e ( m l s ) p r i n c i p l ef o rc o n s t r u c t i n g t h ee f g s h a p ef u n c t i o na r ei n t r o d u c e dt ot h ee mc o m m u n i t y s o m ee x a m p l e ss h o w e f gm e t h o dh a sm a n y a d v a n t a g e so v e rf e m m e t h o d h o w e v e r ，t h ee f gs h a p ef u n c t i o nd o s en o ts a t i s f yt h ek r o n e c k e rr e l a t i o n ，t h ee s s e n t i a l b o u n d a r yc o n d i t i o n sc a nn o tb ee n f o r c e dd i r e c t l y b u ti t c a nb ee n f o r c e db yan e wm e t h o d e f g f e mc o u p l i n gm e t h o dw h i c ho n l yap a r to ft h ea n a l y z e dd o m a i ni sm o d e l e db yu s i n g e fg a n dt h er e s to ft h ea n a l y z e dd o m a i ni sd i s c r e t i z e db yu s i n gf e m i nt h ep a p e ran e w a l g o r i t h mt oc o u p l ee f g w i t hf e mi sp r o p o s e d t h er e s u l t ss h o wt h a tt h eb o u n d a r yc o n d i t i o n s a r ea p p l i e da sw e l la sf e m i nt h ep a p e r ，m e s h l e s sn o d e sa r ed i s t r i b u t e db yb a c k g r o u n dm e s hm e t h o d a n da n e wt h e o r y d o m a i nt o p o l o g yw h i c he f gn o d e sc a nb ei n s e r t i n gr a p i d l ya n d d o m a i ni n t e g r a lc a nb ed o n ee a s i l yi sd i s c u s s e d u s i n ge f g m t os o l v ee l e c t r o m a g n e t i cp r o b l e m ，t h er e s u l t e dc o e f f i c i e n tm a t r i xi s s p a r s e ，s y m m e t r i cb u tn o n p o s i t i v e s oc a n tb es o l v e db yi c c g t h ep a p e rp r e s e n t s an e wn u m e r i c a lm e t h o d l s q rm e t h o dt os o l v e a n dn u m e r i c a le x a m p l e ss h o w t h a tt h em e t h o dg i v e ss a t i s f y i n ga c c u r a c y i naw o r d ，t h ea u t h o rm a d eas t u d yo na p p l y i n gm l mi ne l e c t r o m a g n e t i c c o m p u t a t i o n a n dp r o v i d eo t h e rr e s e a r c h e r sn e wr e s e a r c h i n gt o p i c s k e y w o r d ：e mc o m p u t a t i o n ；e f g m ；e f g f e mc o u p l i n gm e t h o d ；d o m a i nt o p o l o g y ； l s q r 华南理工大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所 取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任 何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡 献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的 法律后果由本人承担。 作者签名：t - 4 - - 越j 钆日期：嘲绰年，月? f 日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意 学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文 被查阅和借阅。本人授权华南理工大学可以将本学位论文的全部或部分内 容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存 和汇编本学位论文。 保密口，在年解密后适用本授权书。 本学位论文属于 不保密曲。 ( 请在以上相应方框内打“d ”) 作者签名：林越凡 导师签名： 第一章绪论 1 1 课题背景和意义 第一章绪论 本课题是受国家自然科学基金委员会赞助的( 批准号5 0 1 0 7 0 0 3 ) ；为广东省 自然科学基金项目( 项目编号9 6 0 2 1 6 ) 的继续研究。 有限元法( f e m ) 是现代重要工业领域中广泛使用、取代具体试验和测试、 以了解物理场分布的主要仿真工具，它可以极大缩短新产品的设计开发时间、节 约设计费用，对近3 0 年来许多技术领域的指数级发展起重要的作用。但是f e m 的技术核心、编程和使用难点以及计算的瓶颈都在于网格的生成和优化( 需要大 量的人工干预、耗费大量计算时间) ，它降低了f e d 使用的方便性和效率，增加 了推广应用的难度；尤其对应用f e m 进行结构优化、动态仿真的场合，由于网格 必须反复生成，需耗费大量的人力和机时。总之，f e m 还存在以下几个问题： ( 1 ) 对于某些问题，f e m 要求对区域反复剖分，这成为计算的瓶颈； ( 2 ) f e m 仅提供c o 的近似解答，当使用标量位或矢量位求解电磁场问题 时，所获得的场的结果沿单元的边界是不连续的； ( 3 ) 网格单元的形状对f e m 的收敛存在严重的影响。 因此，机械材料研究领域的研究者在近二十年中开发了无网格法( m e s h l e s s m e t h o d ，简称m l m ) 。m l m 是继f e m 之后新一代的数值计算技术。它利用覆盖 求解区域的散布点取代f e m 中的网格单元，不仅散布点很容易布置，而且在计 算中不存在保持单元边界同区域边界( 媒介交界) 一致的问题、反复计算不存在 网格再生和优化的问题、计算结果不存在受网格形状影响的问题，解的连续性好、 精度高、收敛快。因此特别适用于形状优化、材料变形、裂纹成长预测等动态问 题以及具有奇点的问题。 在近几年，m l m 已经成为继f e m 后新的研究热点，但研究主要集中在机 械材料领域。在电磁场数值计算领域、特别是优化问题( 如电磁设备结构的优化、 电磁屏蔽装置的优化) 、奇点问题( 如导体尖端电场、电磁兼容的偶极子近端场 问题) 、磁路( 电路) 部件运动问题( 如电磁轨道炮、旋转电机) ，也存在用m l m 取代f e m 、减少计算的人工干预强度、缩短计算时间、提高计算精度的要求。 1 2m i _ m 简介 ( 一) 实现过程 根据e o f i a t e 的定义：任何方法，如果其近似解能够严格按照节点来构建， 华南理工大学硕士学位论文 就是m l m 。即该法基于节点的信息，不需要形成单元或网格。文中用以下分步 骤的操作来解释m l m ，图1 1 让读者获得了赢观印象。 臼日魑日- ，7 i 。、 圈日四日圜 图1 1m l m 实现图解 f i g 1 1m l m 1 ) 在求解区间分布一系列的节点。布点的方式可以均匀，也可不均匀。 2 ) 为每个节点分配一个影响域( 支持域) 。影响域以节点为中心，通常是圆 形或矩形，不同节点的影响域大小可以相同，也可以不相同。 3 ) 在每个节点影响域上，分配一个权函数。它是高阶连续的正凸函数，关 于节点对称，在影响域的边界上及外面等于零。 4 ) 为每个节点分配一个未知参数( 节点自由度) ，并按照某种规则，利用前 分配的权函数，构建形函数，将近似解表达为： 卫 u 5 ( x ) = ：矿，( 工) u ， ( 1 1 ) 五 这里，u ，为节点j 处的自由度参数；n 为区域内分配的总节点数；c k x ) 为节 点j 的形函数。由于形函数构建方法的不同，形成了各种不同的近似方案。 5 ) 对偏微分方程进行离散，得到方程组。常见的有以下几个离散的方案： a ) 点配置法：要求方程在一组点上成立，从而得到相应的离散形式。该离 散方法直接，求解过程很快，但稳定性和收敛速度不理想。 b ) 伽辽金法，该方法基于变分原理。其离散方程的推导和计算中涉及到积分， 通常有3 种方法： ( i ) 单元积分 这是原始的e f g m 在完成区域内积分时采用的，将区域划分成积分单元， 在每个积分单元内采用高阶高斯积分。对于某些积分单元内部存在不连续界面或 区域的边界看似比较粗糙，其实对结果的影响是很小的。 ( i i ) 利用有限元网格作为背景积分网格 2 第一章绪论 这种方案对于e f g m f e m 耦合法是有很大的意义。当划分有限元近似区 域和无网格近似区域后，背景单元在有限元区域作为单元划分，在无网格区域可 作为积分计算单元。 第一种和第二种积分方案最大的不足是没有彻底抛弃单元，但是背景网格 不影响无网格方法的本质 ( i i i ) 点积分 b e i s s e ls ，b e l y t s c h k o t 介绍了e f g m 的一种点积分形式。其作法是：在能量 泛函中增加一稳定项( 该稳定项中含有一个参数一一稳定系数) ，得到相应的控 制方程，从而将区域内积分转化为节点上的积分，该法是可行的。但计算过程复 杂，而且还要选择稳定系数，稳定系数的选择对计算结果有较大影响。文【3 】根 据b e l y t s c h k ot ，k r o n g o u zy , o r g a nd ，e t a 等人的建议：在小积分单元内采用低阶 高斯积分，比在大单元内采用高阶高斯积分得到的解答好，故而提出当积分单元 划分得足够小时，可将每个积分单元内积分用其中心点上的积分代替。 该法积分速度快，而且完全抛弃了网格，但类似与点配置法，可能不稳定。 而且在实现过程中，取多少个积分点才能达到最优? 以及对计算精度是否有影响 尚需探讨。 6 ) 处理边界及交界条件。 这是无网格方法的另一个难点。 ( a ) 边界条件的施加 由于求解域的每点被许多节点的影响域覆盖，因此节点的近似解不仅取决于 本节点的未知数，还同影响域覆盖本节点的那些节点的未知数有关。这一特性可 以解释为形函数不满足k r o n e c k e r 条件： f 1t 一7 仍( - ) 如2 j：- - a j 1 2 也就是说，在x i 处，由式( 1 1 ) 得到的“6 ( 而) u i ，因此第一类边界条件 无法通过类似f e m 的方法直接施加，需另想办法。这类方法很多，有拉格朗日 乘子法，罚函数法，修正的变分法，m l s f e m 耦合法等。 其中通过拉格朗日乘子引入基本边界条件，精确，原理简单易懂，它的最大 的缺点是引入了新的未知量，并且使离散方程的系数矩阵不再具有正定、带状的 特点，适用于规模小的问题。修正的变分原理是将拉格朗日乘子用相应的物理量 代替，这样可以避免由乘子未知量产生的不良影响，同时能保持矩阵带状的特点， 减少了计算量，但是精度低，使用不方便。罚函数法近似于插值法，可达到与 f e m 一样的精度，可用于不规则和规则的求解域，稳定性好，但限制条件很多， 如罚函数最好选择样调函数，且影响域要足够小等，使用起来比较麻烦。而 华南理工大学硕士学位论文 e f g m f e m 耦合法是具有很大的优点，但这种方法在交界区混合形式的形函数 复杂，而且近似函数的导数可能不连续。 ( b ) 交界条件的施加 一般，在不同材料的交界面上，解的法向导数不一定连续，但m l m 提供 c 4 ( d 1 ) 的连续性( 即如果权函数是c 4 连续的，那么形函数和解在全局也是 c 4 连续的，因为通常权函数只可能在影响域的边界处发生不连续，但实际上， 在该处权函数及其任意阶导数的值通常为0 ，所以影响不到形函数及解的连续 性) ，使得法向导数总是连续的，用拉格朗日法和跳跃函数法。 ( 二) m l m 的实现方案 根据不同的形函数构造方法、不同的离散方法，就可以得到不同的m l m 。 具有代表性的有s p h 、d e m 、m l s 、p k m 、e f g m 等。可归纳为3 种基本的实 现方案： 1 ) 函数近似方案 最初主要用于s p h ( 光滑粒子法) ，它对函数利用核函数进行近似。s p h 最 初用于无边界的天体问题。 2 ) m l s 近似 d e m ( 扩散单元法) 是将移动最小二乘方法( m l s m ) 和伽辽金法( g a l e r k i n ) 相结合的方法。 e f g m ( 无网格伽辽金法) ：改进完善了d e m ，这类方法比s p h 计算费用高， 但具有较好的协调性及稳定性。 3 ) 单元分解法 单元分解法是由d u a r t e 和o d e n 等人发展起来的。对于求解区域，单元分解 法用一些相互交叉的子域来覆盖，每个子域都与一个函数m ，( 工) 相联系，函数仅 在子域内非零，并且满足单元分解条件： e 中，( 工) = 1 ( 1 3 ) ， 采用这种近似法的m l m 有：h p c l o u d s ，p u f e m 1 3m l m 发展 m l m 的思想源于8 0 年代。早在1 9 7 5 年，p e r r o n e 和k a o 开发了广义差分 法( g f d ) ，它无需明确的网格却可以处理任意区域。1 9 7 7 年，l u c y 首先提出了 平滑粒子水平法( s p h ) ，用于模拟星爆、宇尘等宇宙无界现象。1 9 8 0 年，l i s z k a 和o r k i s z 提出了单元粒子法( p a r t i c l e i nc e l l ) 。 4 第一章绪论 直到1 9 9 2 年，b n a y r o l e se t 。a l 在求解边值问题的偏微分方程数值解时，首 先用m l s 方法构造形函数。这样形函数的形成及区域积分的实现都可以脱离单 元的概念。bn a y r o l e s 等称这个方法为d e m ( d i f f u s i o ne l e m e n tm e t h o d ) 。1 9 9 4 年t b e l y t s c h k o 等对此做了进一步的改进，改进后的d e m 称为e f g m ( e l e m e n t f r e eg a l e r k i n m e t h o d ) 。1 9 9 5 年，l i u ，j u n 和z h a n g 提出了用修正函数改进 s p h 的方法，与此同时，d u r a t e 和o d e n 提出了h p 云法( h d c l o u d ) ，b a b u s k a 与m e l e n k 对m l m 方法提出了新的阐释，即m l m 方法的“1 ”划分特性。从此 之后，很多的研究者和研究机构开始了对m l m 的研究。 而第一次把m l m 应用在e m 计算中是在1 9 9 3 年，由法国的y m a r 6 c h a l 等 人完成的，它应用了d e m 对f e m 获得的位的近似解进行后处理，得到在整个 求解区域内光滑，高精度的场值。这是应用普通的f e m 插值是不可能做到的， 因为f e m 无法保持解的导数沿单元边界连续。在1 9 9 8 年，他讨论某些m l m 在 e m 计算中的应用问题，同时，v c i n g o s k i 等人介绍了使用m l m 计算电压互感 器电场的研究，a l u r n ，n r 介绍了应用点配置的m l m 于微电子及微机电子设备 的研究；1 9 9 9 年，s a v i a n a 等人介绍了使用m l m p k 方法计算2 d 磁稳态问题 的研究；c h a r a u l t 等人分析了m l m 应用于e m 计算时，边界和交界条件的处 理。 2 0 0 1 年，t a n gy i m i n ，x i ex u e p i n g 将有限元法和小波法结合，提出了小 波一有限元法【4 2 】，它结合了两者的长处，用有限元法施加边界条件，用小波法来 求解场梯度变化剧烈的区域；a l b e r t oc a r p i n t e r i 等则提出了扩展拉氏e f g 法 ( a l e f m ，a u g m e n t e dl a g r a n g r i a n e l e m e n t f r e em e t h o d ) 4 3 】，它通过引入大量 的拉格朗日乘子，却不增加变量的个数，不仅解决了一般无网格法在旌加边界条 件上的困难，而且使边界条件的施加精度得到提高；z h a n gj m 和y a o z h 则提出 了r h b n m ( r e g u l a rh y b r i db o u n d a r yn o d em e t h o d ) 方法【4 4 】，它是基于修正的 变分原理和m l s 原理构建的，同样解决了无网格法边界难于施加的问题。 2 0 0 2 年，g r l i u 提出了p a m 方法( p o i n ta s s e m b l ym e t h o d ) 【4 5 】，它在求解 域的边界上布置一系列的离散点，在域内配置一系列分散节点，这些节点的影响 域由它周围点组成的一系列三角形组合而成。形函数类似于传统f e m ，因此满 足k r o n e c k e r 条件，边界条件可以像f e m 一样直接施加，且计算精度高，无需 背景网格。 2 0 0 3 年，t y n g 等提出了h m d o r m 方法( h y b r i dm e s h l e s s d i f f e r e n t i a l o r d e rr e d u c t i o nm e t h o d ) 【4 3 】，通过把真正的无网格配置技术和核重构近似法相 结合，解决了无网格法边界条件难于施加的问题。 以上这些m l m 法的改进，都是在机械材料领域里实现的，虽然还很少在电 磁计算中得到应用，但这些改进，无疑给e m 计算提供新的工具，为e m 计算方 5 兰要里三查兰堡主堂些笙苎 法的研究提供了新的课题和机会。 1 4 主要研究内容及安排 本文对无网格方法进行了比较全面的研究，包括e f g m 、e f g f e m 耦合法 以及无网格节点的配置问题等等。并且大量列举了静电场的例子，用实例验证该 法的可行性和优势。为电磁场数值计算提供了比f e m 更为实用、方便的工具， 也为电磁设备( 如电磁屏蔽装置) 结构设计、电磁场的奇点，以及电路( 磁路) 部件运动的问题提供了灵活、快速的数值仿真工具与实践经验。论文内容安排如 下： 第二章详细介绍e f g m 的基本思想，以及逐步推导了该法在多媒介电磁场问 题的实现过程。具体讨论e f g m 求解单媒介电磁场问题，以及程序的实现过程； 最后通过一维和二维的实例验证说明该法的可行性和优越性。 第三章介绍了基扩充无网格方法的基本原理以及在电磁场应用，并用实例验 证了该法。 第四章简单介绍了有限元法的基本理论以及在多媒介电磁场问题中的应用， 接着比较了有限元法和无网格方法的异同点，给予了数值计算方法新的启迪。 第五章首先介绍v l a t k oc i n g o s k i 提出的传统耦合方法，并提出了一种新的耦 合法则，它以及在电磁场数值计算中的理论和程序的实现，最后通过实例验证该 法的优势。 第六章实现了用背景网格法来快速配置无网格节点，并且提出了一个新的理 论一一区域拓扑结构化思想，既能快速插入无网格节点，又能实现快速积分。 第七章通过分析无网格法形成的大型稀疏矩阵的特点，其系数矩阵具有稀 疏、对称但不适定，不能用常规的i c c g 法求解。提出采用l s q r 算法来计算， 不但提高计算的速度和进度，同时节省了内存空间。 6 第二章e f g m 第二章e f g m 无网格伽辽金法( e l e m e n t f r e eg a r l e r k i nm e t h o d ，简记e f g m ) 是最典型的无 网格方法之一。它采用移动最小二乘法( m o v i n gl e a s ts q u a r e ，简记为m l s ) 构造 形函数，用伽辽金准则获得离散偏微分方程。同其它方法相比，具有数值稳定、 后处理方便、精度高、收敛快的特点，因而得到了大量研究和应用。凡研究、应 用e f g m 的文献，都要用一定篇幅来介绍移动最小二乘法原理。 2 1e f g m 2 1 1 m l s 基本原理 在域q 中，在节点x f f l = 1 , 2 ，n ) 给定数据“；，则由数据“；可拟合出函数 摊( 工) = p f ( x ) n ，( 工) = p 7 ( m ( 工) ( 2 一1 ) 持l 式中p ( 工) 是完备多项式基，例如： 线性基： p 7 ( ) = l ，工】 p 7 ( 删= l y ) 平方基： ；：兰：叁耋；2 ：，捌，y ：) n ( z ) 是一个依赖位置x 的全局系数，是未知的， 口( 上) = 【( 工) ，嘞( 工) ，- 一，c k ( x ) 】7 1 d 2 一d 1 d 2 d ( 2 2 ) ( 2 3 ) ( 2 4 ) 12x34 图2 一l 覆盖点j 的影响域 f i g 2 1o v e r l a p p i n gd o m a i n s o fi n f l u e n c ea n d l o c a ln o d en u m b e r i n ga tp o i n tx 为确定任意点工处a ( 工) 的数值，将“( 并) 在工周围局部近似成以下形式( 即固 定a ( 工) ) ： “( 膏) ；“( 孑) ，= p j ( - 2 ) a i ( x ) = p 7 ( i ) a ( z ) ( 2 - - 5 ) 华南理工大学硕士学位论文 其中h 表示近似的适用范围，由覆盖点膏的那些影响域构成，如图2 一i 所 示。上式在i = 工处是精确成立的。 借鉴最小二乘法曲线拟合的思想以确定a ( x ) ：在“一( i ) ，近似中，在那些影响 域覆盖点x 的节点研处，近似值“6 ( x ，) 。和样本值( 自由度值) 蜥存在偏差，使这 些偏差的平方和达到最小，可以确定a c x ) ；考虑到影响域覆盖x 的那些节点工， 中，离工较近的，应该对a ( 工) 的贡献较大，离工较远的，应该对a ( 上) 的贡献较小， 故在局部范围内采用加权误差平方和最小的策略。具体操作如下； j ( x ) = w ，( x x 。) p ( x ，x 。) 一“：】2 ：兰m ( x - - x 1 ) p r ( x 。) a ( x ) 一“；】z ( 2 6 ) = p a ( x ) 一u 4 r w ( x ) t p a ( x ) 一u 】 其中w i ( x x i ) 是节点i 的权函数，它只在以x i 为中心的局部区域( 称为x i 的 影响域) 内非零。 式中 p l ( 而) p m ( j 1 ) p 7 ；i ； i i ， 【p l ( j ) p 。( 确) j w l ( x x 1 ) w ( x ) ：l o 【 ； ( 2 7 a ) ( 2 7 b ) ( 2 7 c ) i ，( 上) 中的求和虽然是对所有节点x ，进行的，但只有影响域覆盖到点x 的那些 工，其w ，“x ，) 才不为零，故，( 工) 只由x 点附近的一些节点决定。 令塑盟= 0 ，可得： 3 a ( x ) 享墼：p t w ( x ) 【p a ( x ) - - u $ 】- o ( 2 8 ) d a ( x l 等a ( x ) = p w ( x ) p 7 】- 1 p w ( x ) u 记： a ( 工) = a 一1 ( 工) b ( j ) u ( 2 - - 9 ) 其中： a ( 工) = p w ( x ) p 7 第二章e f g m 最后，令i = x ，得到 其中形状函数为： b ( x ) = p w ( x ) u ( x ) = u h ( i = x ) 。 = n ( 工) q ( 石) = p r ( x ) a _ 1 ( 善) p 【w l ( x x o u i w ( 工一h ) “ 1 = o a x ) w a x x j ) u ， i = 1 - - z o a x ) u ， ( 2 1 1 ) ( 2 1 2 ) 嘭( x ) = 珥( z ) h ( 工一) ：( x - x , ) 缶7 ( ) a 一1 ( 工) p ) ， ( 2 1 3 ) ：w ，( x - - x t ) 兰b ( 工) a 一1 ( 算) p l ， j = l 这种确定形状函数的方法称为移动最小二乘法。“移动”的意思是指：系数 a ( x ) 同位置有关，即为位置的函数，将点x 移动到某处，以上方法就可以确定该 处的a 佃) ；而在传统的最小二乘法中，系数则固定与位置无关。 另外，由于计算需要，还要确定形函数的导数，文中分别就一维和二维情况 讨论形函数的求导问题。 ( 一) 一维 取： p 7 ( 工) = 【1 ，x 】 ( 2 1 4 ) 则： a ( j ) = w ( x x ，) p ( x z ) p 7 ( _ ) ：石一 ， 二 + “x 一工：， 乏 十+ w c x 一k ， 2 1 5 a 又 口( 工) = 【w ( x x o p ( x o ，w ( x 一屯) p ( ) ，“ = w , - i x x 州一柑 u 刀胍卜k 飞 肛 妒 华南理工大学硕士学位论文 办。= ( p 7 a 。1 马) ， = p r a b f + p 7 ( a 一1 、，b f + p r a b ， 其中： 哆，= i d w ( x 一而) p ( _ ) a 一1 j = - a 一1 a z a 一1 这里： a x ( 石) = w ，( 工一x t ) p ( x 1 ) p ( 而) ：警1 = 1 。一， 二群x - 1 + d 出w 。一屯， 也x ；1 j 坐d x 。一矗，1 2 1 9 ( 二) 二维 取： p 7 ( x ) = 1 ，x ，y 】 ( 2 2 0 ) 则： a ( x ) = w ( x x ，) p ( x ) p 7 ( x j = w ( x x 1 ) 1 x i 砰 咒彬 曰( x ) = 【w ( x x 1 ) p ( x 1 ) ，w ( x - x 2 ) p ( x 2 ) ，w ( x - x 。) p ( x 。) 】 = w c x - x ，f ，w x - x ：， ，w c x - x 。， 1 y n x 。2x n y 。l y nx ；y ，y 2 。1 ( 2 2 1 ) 其中x = k ) ，】和x 。= x a ，y 。】 一维时形函数的导数比较容易求得，但是二维就很费事了，本文采用 b e l y t s c h k o e ta l 提出的方法来求解形函数的导数，该法对a 用l u 分解。形函 数( 2 - - 1 3 ) 可重写为： 办= p r ( x ) a 。( x ) 岛( x ) = f ( x ) b f ( x ) 则： a ( x ) y ( x ) = p ( x ) a 采用l u 分解式p ( x ) a ( x ) = l ( x ) u ( x ) ，上式两边同乘以尸( x ) ，得 1 0 k x “ + 、i1 m 秽建屯w 屯m l k x 以 + 1，j m u 井 第二章e f g m j 玳u 訾1 2j 匕：訾 jl岛u。1：1：75：+要ua：2y：2-：-笼u1：37：3：乃73，)+k(u2m27。2托+u2373)121j r u 2 3 7 3 + u 2 3 7 3 ，+ 坞，为 = 圣 = l ( “l l + m + “1 2 托+ 3 l = ic 2l i 岛l ( l4 x + q 2 虼+ “1 3 乃) + f 3 2 m 2 2 托) + 坞3 为jl c 3 j 这样就把州x ) 求出来了，接着就可以求形函数的导数了。 艮= 【o 1 0 】 中，( x ) = 广，( x ) 易( x ) + 广( x ) 目，。( x ) 2 1 2 权函数 权函数的选择很重要，要遵循以下原则”1 ： 1 1 具有紧凑的定义域( c o m p a c ts u p p o r t ) 2 ) 非负，且在指定范围以外为0 ( 2 2 5 ) f 2 2 6 ) ( 2 2 7 ) r 2 2 8 ) ( 2 2 9 ) r 2 3 0 ) ( 2 3 1 ) f 2 3 2 ) 1j q 吃q ，l 1 1 1j m 托乃 。l 术 1，l，j u ” u u u 华南理工大学硕士学位论文 3 ) 在x i 的近处大，远处小 4 ) 至少一阶和二阶偏导数存在( 即高阶连续性) 权函数的选取没有理论上的具体规则，带有某种任意性，可选用指数函数 锥形函数，三角函数，样条函数等。如可选用三次祥条函数如图( 2 2 ) ： w ( x 一再) = w ( r ) = 三一4 r 一4 ，3 3 兰一4 r + 4 ，一兰， 33 0 r 三 三 r 1 其中r = 乡乞，其中d ，= i i x - x , i i ，是z 与x ，的距离；d m 是x - 的影响半径 ( 2 3 4 ) 其中，一般取2 d 一4 。 c ，：m a x 忙一x 川( 2 - - 3 5 ) ，es t 式中s f 是x l 的邻近节点组成的最小点集( 能保证形函数的表达式中矩阵可求 逆的最小节点数目) ，c ，表示s ，中节点到x ，的最大距离。 对( 2 3 3 ) 式求导，得： 地：业生： 出咖出 ￥ ( 一8 r + 1 2 r 2 ) s i g n ( x x j ) ( 4 + 8 r 一4 r 2 ) s i g n ( x - x ，) 0 图2 2 三次样条函数 f i g 2 2c u b i cs p l i n ew e i g h t f u n c t i o n 1 2 ( 2 3 6 ) 弓扣 茎三童婴翌 2 2 单媒介静电场边值问题求解 2 2 1 离散格式 本节采用另一种方法一一修正的变分法来推导其离散格式- ”。 考虑以下二维电磁场问题： 邮等十哮一 f d ：妒= l ：娑= 一吼 d ，l 上式相应的弱条件的能量泛函为： 脚) = 皿 譬【( 誓2 + ( 等2 卜即 撕一( f 。卢警舾一j ：r n q o 油) ( 2 - - 3 8 ) 用拉格朗日乘子法施加边界条件，则对应的修正弱条件的能量泛函为： f ( 妒 五) = f ( 咖+ 【五( 妒一) 豳 ( 2 3 9 ) 这里五是拉格朗日乘子，其表达式为： a ( j ) i r = n ，( s ) 丑 ( 2 4 0 ) ( 其中n d s ) 是一类边界上的拉格朗日插值基函数。丑是待求系数，s 是边界弧长) 求( 2 - 3 9 ) 式的最小值，得到修正的弱条件变分方程： 万f f y 伊，卫) = j 瞧f f ( g p ) ，织+ ( 巧) ，妒，卜a g p d x d y ：+ 。j ：v q 0 8 q j d s - l 。卢鼢出+ 。栅 ( ：叫， 万，( 妒，五) 2f 。觑7 ( 妒一# o o ) , i s = o 其中 将( 2 - - 4 0 ) 式和( 2 1 1 ) 式等代入( 2 4 1 ) ，并且令却= m 。( x ) ，觑= n j ，可得： 臣孑胁 f ) ( 2 - - 4 2 ) 卜 华南理工大学硕士学位论文 f ，j 2 j j p t ( , i , j ) ，( m ，) ，+ ( m ，) y ( m ，) y a x a y ，为方阵k 的元素； g = lm 。n j d s ，为矩阵g 的元素； 工：见砷，a x a y 一工。吼m 。d s ，为列向量f 的元素； 2 4 3 t i = ln ；“o d s ，为列向量t 的元素； 当然( 2 - 4 3 ) 式也可以用伽辽金法推导得到。 推出e f g m 的离散格式后，应用高斯积分等数值积分方法，计算以上各系数 矩阵，再求解方程组( 2 4 2 ) ，可得数值解u ，由u 按m l s 拟合出的便是问题的 数值解。 2 2 2 程序简介 用v b 语言，结合m a t l a b 公司提供的矩阵计算v b 插件m a t r i x v b ，进行了 程序实现。其程序流程如下： 1 。产生m l m 的节点； 2 。由边界信息确定背景网格，进而确定域内高斯点坐标、所属区域； 3 。将边界划分为高斯积分单元线段，并确定边界高斯点坐标； 4 。建立域内高斯点一节点关联矩阵； 5 。建立边界高斯点一节点关联矩阵； 6 。对域内高斯点扫描，计算对于区域积分的系数； 7 。对边界高斯点扫描，计算对于边界积分的系数； 8 。利用m a t r i x v b 的矩阵运算功能求解方程组，获得各节点自由度参数； 9 。对任意点，利用高斯点一节点关联矩阵，快速确定其与节点的关联关系， 再计算其位值、场值。 以1 d 泊松方程的m l m 实现过程为例，简单介绍程序的编写流程( 采用 m a t l a b ) ，虽是l d 问题但很容易推广到其它类似问题。程序首先确定有关常数和参 数，如权函数的参数d 设置为2 。接着布置节点。在域内分布n 个节点( 由这些节 点的电位，可拟合出区域内其他点的电位) 。这些节点可以均匀分布或不均匀分布。 为使程序简化，在【0 ，1 】内均匀布置“个节点，这样各节点的c 。和d l i l i 都相等。 而后进行区域积分。首先把求解域分成几个积分单元( c e l l ) ，在每个积分单 元上应用高斯积分法求解方程中的各积分系数。对于每一个积分单元，先确定高 斯点。为简单起见，程序中取了1 0 个积分单元，各个积分单元内只取一个高斯 点( 两相邻节点的中点) ，则每个积分单元的雅克比系数j a 都等于o 5 ，存放在矩 阵g g 中。然后对每个高斯点，找出影响域覆盖到它的节点，并计算该高斯点与 这些节点的距离d i f 。接着计算节点在该高斯点处的权函数及其导数。注意，对 于影响域不覆盖该高斯点的节点，相应的权函数及导数为零。再计算该高斯点处 1 4 第二章e f g m 的形函数及其导数，具体过程为：先计算b = p 1 w ( x ) ，存放在矩阵b 里；再计算 a = p w ( x ) p ( x ) 及其导数，分别存放在矩阵a 和a i n v 里；而后按巾= p r a 1 ( x ) b ( x ) 计算权函数及其导数，结果分别存放在p h i 和d p h i 里，其中p = f 1 ，x ，其导 数为 0 ，1 。形函数计算中牵涉到矩阵的求逆，而m a t l a b 语言对矩阵求逆有很 多方法，如l u 分解法。程序中直接用函数i n v 0 来求矩阵的逆阵。然后用式( 2 - - 4 3 ) 计算k ，f 中的第一项。以上计算应遍历完所有的高斯点和节点。 紧接着进行边界积分。其过程与区域积分类似，主要是计算n ，g 以及f 中 的第二项。在此不再赘述。 经过上述过程可得到离散方程中k ，n ，f ，g 值。最后求解这个离散方程，得 到n 个节点的自由度参数u 。，由式( 2 4 2 ) 可得各节点的电位。 2 2 ，3 算例 2 2 ，3 1 一维算例 考虑一个无穷大平行板电容器的静电场i 司题。设该电容器的两极板之i 刖充有 密度为p = e 的自由电荷，且两极板都接在+ 0 5 v 电源端，极板问的距离为2 ， 如图2 3 所示。由于电容器的激励和几何形状都关于y 轴对称，故而只要求解 一半的区域，另一半的电场分布可以从对称关系得到。为简化问题，不考虑实际 的物理单位。其精确解为：睇( z ) = 1 一x 。这是一个一维静电场问题，可描述如 下： 厂 v 2 “= 一1 x ( o ，1 )( a ) jr d ：u l = 0 5 ( b ) ( 2 4 4 ) l k ：乳= 。 0 5 v0 5 v l r lo + l o u ：o d x 图2 3 两极板电容器示意图 f i 9 2 3c a p a c i t o r sf i e l d x 华南理工大学硕士学位论文 ( 一) e f g m 计算结果与精确解的比较( 采用1 1 个均匀分布节点，1 0 个积分单元) 栩确l i 和e f g m 解的比较 图2 3 精确解与数值解的比较 f i g 2 3r e s u l t sc o m p a r i s o nb e t w e e ne f ga n de x a c to n e 可以看出，用e f g m 计算的结