（电工理论与新技术专业论文）棱边有限元法的理论研究及其在工程涡流问题中的应用.pdf

声明 本人郑重声明：此处所提交的硕士学位论文棱边有限元法的理论研究 及其在工程涡流问题中的应用，是本人在华北电力大学攻读硕士学位期间，在 导师指导下进行的研究工作和取得的研究成果。据本人所知，除了文中特别加以标 注和致谢之处外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为 获得华北电力大学或其他教育机构的学位或证书而使用过的材料。与我一同工作的 同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 签名：王恙帆日期：剑 关于学位论文使用授权的说明 本人完全了解华北电力大学有关保留、使用学位论文的规定，即：学校有权保管、并向有关 部门送交学位论文的原件与复印件：学校可以采用影印、缩印或其它复制手段复制并保存学 位论文；学校可允许学位论文被查阅或借阅：学校可以学术交流为目的，复制赠送利交换学 位论文；同意学校可以用不同方式在不同媒体上发表、传播学位论文的全部或部分内容。 ( 涉密的学位论文在解密后遵守此规定) 作者签名 豆志、日日 导师签名：摊 日期：竺三! ! ：!日期：型： 华北电力大学硕士学位论文 1 1 计算电磁学发展简介1 】 第一章引言 计算电磁学( c o m p u t a t i o n a le l e c t r o m a g n e “c s ) 是一门新兴的边缘交叉学科一， 它以电磁场理论为基础，以高性能计算技术为手段，运用计算数学提供的各种方法， 解决复杂电磁场理论和工程问题，电磁场与微波技术学科中一个十分活跃的研究领 论。 计算电磁学的研究与发展可以分为电磁场分析( 正问题) 、逆问题求解( 含优化 设计问题) 和电磁场工程三个部分，它们相互衔接，又相互融合、相互促进。近几年 来，电磁场工程在以电磁能量或信息的传输、转换过程为核心的强电与弱电领域中 显示了重要作用。电磁场工程面对的是复杂的大系统工程问题，其中常包含电磁场 及相关物理场在内的瞬态耦合问题、优化设计问题和逆问题，通常还含有非线性 数十年来，随着计算机技术的飞速发展，计算电磁学领域已经取得了重大的科 学技术进步。其中：有限元法( f e m ) 的二维、三维解已经有了很大发展，包括对稳态、 时变场问题和非线性问题、运动媒质问题的处理，对规范问题的正确理解等等。用 有限元法解决工程问题的论文所占的比例最大。目前，三维涡流场分析仍然是最受 重视的问题之一。边界元法( b e m ) 可以用来分析二维、三维问题，但边界元与有限元 相比较，那一种方法更有前途，仍没有定论。用棱边元( e d g ee l e m e n t ) 方式构造f e m 和 b e m 的基函数，是对传统节点元的革新，对于描述场的变化和连续性提供了有效的物 理框架。 高性能计算技术的发展是2 0 世纪后半叶最重要的科技进步之一，计算电磁学 是其中的重要分支。它不仅是现代电工、电子和信息技术的理论基础，而且也成为 军事、生态、医疗、天文、地质等众多领域新技术理论的生长点。与此同时，现代 科学技术各学科之间的相互交叉、渗透，工程技术日趋集成化的特点，正在进一步推 动计算电磁学的发展。 1 2 涡流问题数值计算发展概述【2 , 3 , 4 目前，二维场及三维静态场的分析已基本成熟，而对于三维涡流场虽然人们 对它进行了大量的研究，但远未达到成熟的阶段。三维涡流场的分析较之静态电磁 场的分析复杂，除了电磁场随时间变化外，还涉及多连通域、材料的非线性、磁滞 等等因素的影响。因此，对三维涡流求解方法的研究仍然是工程电磁场领域的前沿 课题。 八十年代初，各种涡流问题分析方法相继发表，从方法上可以分为：边界元 1 华北电力大学硕士学位论文 法，有限元法，有限元和边界元混合法等。 边界元法是把积分形式的控制方程在边界上离散得到一个代数方程组，这比有 限元法对同一问题有较小的，但不稀疏的系数矩阵。当边界上的未知数求出之后， 再求场域内任意点的场。 有限元法于1 9 4 3 年被首先提出。它以变分原理为基础，对齐次第二类边界条 件和不同媒质分界面的边界条件不必处理，形成的系数矩阵为稀疏矩阵，因此在解 决非线性，多媒质电磁场问题上有其特有的优势。但是，用有限元法处理开域问题， 有一定的困难，必须人为地添加外边界，使计算场域为闭域空间，若外边界选取太 远，会有计算代价商的问题。 常用的涡流场计算方法，从求解变量上可分为：基于a 有的限元法；基于t 的 有限元法、a 与t 混合的有限元法以及直接以e 或h 为变量的有限元法。此外根据 在涡流区除引入矢量位外是否引入标量位等问题派生出一系列的分析方法，( 1 ) 基 于a 的方法a 一( p a 法，即在全求解场域用矢量位a ，涡流区引入标量位( p ，最 先被提出。该法对多连通域不需另作处理，但未知量较多，求解规模大。若在非涡 流区引入标量磁位掣，取代矢量磁位a ，则为a 一y 法，该法有较少求解变量， 但在涡流区和非涡流区交界面存在修正矢量磁位a 和标量磁位掣耦和的问题，也 增加了问题的复杂性。( 2 ) 基于t 的有限元法如t 一、王，法，导电区域( 包括涡流区) 内选择矢量电位t 和标量磁位t 共同求解，在无电流区域内只有标量磁位、壬， t 一、壬，法计算代价低，数值实现简单，解的稳定性好。但是t v 法不能处理多连 通问题。一个简便的近似处理方法是用电导率接近0 的“导体”材料把形成多连通 域的孔填充，把多连通域转化为单连通域问题。( 3 ) a 与t 混合的有限元法如o b i r o 提出的t 一掣一a + 一v 法和文献 2 提出的a 一( p t 一甲组合分析方法。( 4 ) 直接 采用h 或e 作为求解变量，可降低计算代价、避免微分误差。传统的节点有限元法 求解矢量场h 或e 时，把场矢量的各个分量看成是处处连续的，即使在不同媒质的 交界面也是如此，这和物理实际相矛盾。因此传统的节点有限元法在应用上述方法 处理多媒质涡流问题时将产生困难。当然可以在不同煤质分界面处加大剖分密度使 得求解量在分界面处迅速变化来近似模拟实际情况”3 ，但是这样作显然增大了剖分 难度和计算量。解决这一问题的有效方法是应用棱单元。 1 3 棱边有限元法简述 文献 3 已对棱边有限元法的发展作了详细综述，本节只对其中有代表性的内 容作简要回顾。 1 9 5 7 年，数学家a w h i t n e y 在微分几何引论一书中提出了棱边有限元的概 念，为棱边有限元法奠定了理论基础。1 9 8 0 年，n e d e l e c 首先讨论了混合单元在电 华北电力大学硕士学位论文 磁场和弹性力学中的应用，指出：棱单元由于其切向分量连续而法向分量不强制连 续，可以求解h 或e ，同时推测也可以用于求解矢量磁位a 。1 9 8 2 年法国学者 a b o s s a v i t 提出了以h 为变量求解三维涡流问题，研制了著名的电磁计算软件 “t r i f o uc o d e ”。1 9 8 2 年，a b o s s a v i t 和j c v e r i t e 用有限元和边界元相结合的方法 解开域涡流问题，以h 为状态变量。1 9 8 8 年程志光建立了包括棱单元、节点元和混 合元( 位于导体区和非导体区的分界面处) 的混合有限元模型。1 9 9 0 年任卓翔提出 在棱单元中以e 为状态变量。至此，棱单元法的基本框架从理论和实践上都形成了。 之后的工作集中在解决各方面更细节问题，如：棱单元的插值函数、棱边有限元法 和其它数值方法的耦合、棱单元与节点元方法的比较、规范问题、法向约束问题等 等。 另外需要提及的是虽然棱单元法在一些方面表现出比传统的节点元法优越，被 称为富有革命性的方法，但是棱单元法的发展不论在理论还是在实践上都不及节点 有限元法成熟。一些有争议的问题( 如：同种介质中单元法向连续约束是否需要强 加、棱单元和节点元的比较以及假解问题) 有待解决，商用软件还不如节点元法丰 富。总之，棱单元法还处在一个逐渐发展的阶段，它和节点元方法各有所长，正如 程志光在他的专著【2 l 中所说的那样：“棱单元和节点元是有限元家族中的两个主要成 员，它们的应用卓有成效，两者之间并不存在“谁取代谁”的问题，呈现并行发展 的趋势。而且有趋向表明，两者相结合也许更有前景”。 1 4 本文主要研究工作 第二章论述了棱边有限元法的唯一性问题与树规范，证明了在不加树规范情 况下，b i c g 法可以求棱单元方程组。理论上，若b i c g 法可正常进行迭代，则迭代 次数至多需迭代m 次( 为a 的秩数) 。 第三章对以e 为变量的棱单元法的唯一性和边界条件作了证明并对泛函方程进 行了离散，介绍了电流源区处理方法，推导了基于牛顿一拉夫逊法的非线性迭代 格式。 第四章介绍了程序设计阶段对数据结构方面的改进，实践表明本文提出的数 组与改进十字链表混合的存储技术可以有效地克服了十字链表存储技术的弊端，不 仅节省了约2 5 的内存占用量，而且还“意外”地节省了约2 0 的c p u 时间。 第五章用e 为变量的棱边有限元法分别计算了t e a m 2 1 问题的线性和非线性 模型。 第六章为结论和有待进一步解决的问题。 华北电力火学硕士学位论文 第二章b i c g 法求解奇异方程组的证明 2 1 棱边有限元法的唯一性问题与树规范 在电磁场数值计算中，棱边有限元法作为节点有限元法的重要补充受到了众多 学者的瞩目，它可以克服节点有限元法一些难以克服的困难。在某些方面，棱单元 显示出比相同位函数的节点元优越，但棱单元法同节点元一样也遇到了一些问题， 规范问题就是其中之一。 众所周知，应用节点元以矢量位a 为变量求解三维涡流场时，a 的双旋度方程 和节点元都不能保证a 的散度条件，但是可以用加入罚函数的方法解决。同样，在 应用棱单元以矢量位a 为变量求解三维涡流场时，a 的双旋度方程和棱单元也都不 能保证a 的散度条件，不同的是由于在棱边上a 散度无法定义所以不能用节点元法 加入罚函数的方法加入规范。这一矛盾可以用r a l b a n e s e 建议的定义树一余树的方 式解决。首先根据不同情况，将求解区域内的所有棱边划分为树和余树两个集合， 然后令树棱边上的a 为0 ，只求解余树棱边上的a ，实际上这是一个施加规范的过 程( 关于“树一余树”规范的详细讨论见参考文献【6 】) 。但是“树一余树”规范所 得矩阵的条件数因所选的树而异，故有选最优树的问题。实践证明，在任意场域选 树可以用多种方法，但是选择最优树是困难的【3 。若不加规范方程组的系数矩阵 是奇异复对称阵。有趣的是用通常的一些预处理b i c g 法能得到方程组的解，且已 被用于实际计算中。另外，据文献 7 的推证结果，不论如何选树，所得的余树阵的 条件数都是增大的，导致实际中规范后的矩阵迭代次数增多收敛性变差精度下降。 总的来说，加规范除了能保证方程组的唯一性外并无其它优势。而且在该问题中a 只是中间变量，工程中关心的是磁感应强度b 或电流密度j ，它们是唯一的，a 是 否唯一并不重要。 b i c g 法可以求解非奇异方程组早有证明0 1 。但实践中，该法也已被用于求解 奇异方程组，但理论上还没有关于该法是否可以用于求解奇异方程组证明。本章以 下内容将证明b i c g 法可以求解奇异方程组，而且理论上的迭代次数与树规范后非 奇异阵的迭代次数相同。 2 2 b i c g 法求解奇异方程组的证明 设a 为h e r m i t e 阵，若令b i c g 法的余量为r ，辅助余量为晨，其初值盖。= 占一删 而不是一a s ) + ，则b i c g 法的迭代格式与c g 法迭代格式相同，可把c g 法看作b i c g 法的特例，所以下面的讨论以b i c g 法作为基础( 只在必要时提及c g 法) ，结论同 样适用于c g 法。 华北电力大学硕士学位论文 由于b i c g 法是一种经典的半迭代方法所以本文只证对经典b i c g 法扩展的部 分，其它证明请参考文献 9 。 2 2 1 求解奇异方程组的b i c g 法 记= ) 其中1 表不取转旨，十表示取共轭。若a c ”满足a “= a ，则称a 为 h e r m i t e 阵。记( ，o 为通常复内积，x = i x ，x ：z 。】7 ，y = 砂。y ：y 。 r ，咒c ， i = l ，2 ，n ，则( 置y ) = x ”y = z x ：y , 。 定理1 设矩阵a c ，且爿“= a ( 即月为h e r m i t e 阵) 则向量戈是方程( 1 ) 的解的 充分必要条件为量使式f i x ) = ( a z x ) 一似动一慨习取稳定值。 证明：记 x = k l 而f ，占= b6 2 玩】，x ib ，c ，i = 1 ，2 ，；月 则f ( x ) - - ( x x x ) - ( x , 日伍均 ( 2 一1 ) = l 口，x ，+ l z ；- z b 。x ；+ 一b ，x ； 掣= 毒 喜( 鼢， x 捧y = l 巧 z ， a 慨功+ 慨硎： o x去降i 审薯 = 可 s ， 由式( 2 - 1 ) 、( 2 - 2 ) 、( 2 - 3 ) 式得g r a d f = ( x x b ) 故g r a d f = 0 当且仅当 a x = b 时成立，定理1 证毕。 定理1 中隐含说明了方程躺= b 必须有一个或多个解，但不能无解，也就是说 方程必须是相容的。 定理1 提供了从求解任意n 阶非奇异复阵的c g 法推广到任意求解1 1 阶奇异复 阵的理论基础，因为若a 不是h e r m i t e 阵则可附加一个辅助方程组及其解向量贾和 某指定的右端项画，4 “霄= 占此式与a x = b 联合形成一个c ”c “上的2 n 阶方程 组，其系数矩阵m 是h e r m i t e 阵，即： 匕公 爿0 倦 a ， 对方程( 2 - 4 ) 应用扩展后的b i c g 法既可得到删= 日的解。 定理2 若b i c g 法可正常进行迭代，则迭代次数至多需迭代m 次，m 为a 的秩数 证明：设有一组数屈c ( 1 = 1 ，m ) ，使f r r 一0 ，则有( y 卢i r ，晨，) = o 百莒 。 ( = 1 ，m ) ，利用关于r 的双正交性( 衷，月，) = ( r ，詹，) = 0 i j ，得，( 五，r ，) = 0 ， 又因 ( r ，r ) 0 ，因此卢i = 0 ( i = 1 ，m ) ，因而r i ( i = 1 ，r f l ) 线性无关。又因 的秩为m 且显然r ，a 空间，所以r ( i = 1 ，m ) 是a 空问的一组基。同理第m 次 华北电力火学硕士学位论文 迭代产生的r 与这i n 个向量线性无关，于是r 。= 0 。证毕。 定理2 中“可正常进行迭代”是指不发生下述情况：在m 次迭代内对某个i ( i m ) 可能发生( 乒，爿乒) = 0 使计算，屈时分母为o ，或者发生对某个j ，( 盖，r ，) = 0 ，使 爿。= ，从而迭代中断经典b i c g 法理论上不能对于所有的方程组排除这种潜在 问题需要指出的是经典b i c g 法扩展到奇异方程组并不增加这种潜在的危险( 相对 非奇异的情况) ，奇异阵显然有一个关于a 的1 1 - m 维的零空问( 由r i ( 7 = 1 ，m ) 的线性组合构成的空间记作，的正交补空间记作1 ，维数n m ，若ke n l 则 有a k = 0 ，所以1 为爿的n m 维的零空间) ，但是由迭代公式易知向量 只，甲，r i ，r ”( i = 1 ，1 1 1 ) 都属于爿的空间( 同非奇异的情况) ，故经典b i c g 法扩 展到奇异并不增加这种潜在的危险。另外，若a 为正( 负) 定或半正( 负) 定实 h e r r n i t e 阵时可理论上保证不发生上述情况，此时可用c g 法的迭代格式且有 ( r ，r ；0 ，r j 0 及( 只，爿只) 0 ，只0 。这一点也同非奇异的情况的c g 法。 如果发生上述潜在情况，解决方法同基本的b i c g 法 若a 为非奇异阵，可把原方程变为a ”栅= 4 “曰，显然变形后的方程系数阵为正 定h e r m i t e 阵且与原方程是同解方程组，可用c g 法求解。但对于奇异阵，方程 爿7 a x = a 7 b 与原方程并不同解( 可看出删= 口的解集只是爿”a x = a ”b 解集的 部分) ，用c g 求解此类问题可能得到错解，故应该用b i c g 法求解。 求解奇异阵的线性方程组并非是b i c g 方法的专利，一些迭代方法也可以解，但 并非都可以求解。 以上仅是从理论上对经典b i c g 法进行扩展，实际计算时由于误差的影响迭代 次数同样会超过理论值。 2 2 2 b i c g 法的收敛性 经典b i c g 法的条件数的定义：文献 1 0 指出当a 为非奇异阵时c o n d = 孥生，其 6 t 1 1 j n 中s 、s m i 。分别为最大和最小奇异值。对正规矩阵有s i = n ，故非奇异阵正规阵 j1 c o n d = 鲁爿，其中l k 。小i h ，f 分别为按模最大和最小特征值。文献 8 指出由于b i c g l n 如l 法用于半正定阵时0 特征值对应的特征向量对迭代过程无影响，故可将条件数定义 为c o n d = ； ，中刀。为最小非零特征值。鉴于以上原因可将奇异阵对应的条件数 l m l n a 定义为如下形式：c o n 扛i 2 5 ，其中s 。、s o m 分别为最大奇异值和最小非零 ow l n 1 1l 奇异值。对正规矩阵有s = n ，故非奇异阵正矩规阵c o i l d = 略2 6 ，其中i k 。小 p m * j i 矿m n l 分别为按模最大和最小非零特征值。显然奇异或非奇异方程组的条件数可根 据不同的情况用2 - 5 式或2 6 式统一表达。 华北电力大学硕士学位论文 2 2 3b i c g 法用于电磁场数值计算中的相关问题 2 2 3 1 预处理问题 用有限元法求解电磁场问题最终归结为求解一组稀疏的代数方程组。对于大型 的三维工程问题，有限元离散方程组的阶数可达上万阶，其条件数可能很大，基本 的b i c g 法无法得到收敛解。因此系列预处理b i c g 法被提出并有效的解决了电 磁场数值计算的大型方程组问题。对于用棱边有限元法计算产生的奇异方程组同样 也存在这个问题且用b i c g 法求解时同样可利用预处理的方法。例如文献 3 通过求 解时谐涡流场问题得出了如下结果：当用节点元计算涡流场且方程组非奇异时， p i c c g ( 0 ) 法是首选，当计算规模很大时，可采用p i c c g ( 云= b a b ) 法。当用 棱边有限元法计算涡流场且方程组奇异时，p i c c g ( l = b a b ) 法是首选。 2 2 。3 2 规范问题 以a 法为例，在电磁场的有限元计算中产生奇异阵的原因是未知量x 的散度未 定义，而b i c g 法有解是因为在迭代的过程中a 的散度被弱定义】。具体说：设 g 7 为散度算子的离散形式”】。把g 7 作用到x ( x 为a 法中的解向量) 即g 定 义了弱意义上的x 的散度。因为g 7 p = d 且鲰+ 严甄+ f r ，故有g 斛严g 7 x = g 2 而。这就意味着在每次迭代中j 的散度被弱定义且等于解的初值】。 2 2 4 结论 本文证明了b i c g 法可用于解奇异方程组若b i c g 法可正常进行迭代，迭代次数 不超过矩阵的秩。从本文证明可看出b i c g 法解奇异方程组与解非奇异方程组比较 并无特殊性，因此可把非奇异情况的证明看作奇异情况证踢的特例，b i c c 法求解非 奇异方程组是求解奇异方程组b i c g 法的特例。 由于本文证明了b i c g 法用于解奇异方程组的完备性。对于像以矢量磁位为变量 的棱单元法解三维涡流场的例子，则可用预处理b i c g 法求解无须强加规范条件。 华北电力大学硕十学位论文 3 1 唯一性问题 第三章以e 为变量的棱边有限元法 考虑如图3 1 所示的涡流问题，忽略位移电流的影响并假设其为时谐问题，则 m a x w e l l 方程可写为如下形式： l v h = - ，+ e r e + ，国扰 【v x e = 一j c o , u h ”。曼2 口 s 7 边界条件 虬 l 厅x 日= bo r s 2 方程中e 为电场强度，日为磁场强度，占为磁 通量密度，以为源电流密度，以为涡流密度，为磁 图3 - 1 涡流蒿题求解区域 导率，盯为电导率， 为虚数单位，c o 为角频率。 证：设问题有两组解( e ，日。) ，( e “，h ”) 。 令e = e j e 。h = h j h 。， 则对( e h ) 有： l v x h = j + o e + j c o t e 【v e = ，啤脚 f n e = 0o ? ls 1 i n h = 0 o ns 2 v ( h xe ) = e - v x h h v x e = a e e + i w e e - e + i c o p h h 两侧积分得：h e d s = a + 嬲川e l2 + j c o r i h 7d v 忽略位移电流得：h e d s = p 例2 + j c o u l h t 2 d v 由齐次边界条件可知方程左侧为0 。 由于不为零，全域有h = 0 ，故全域磁场强度唯一。 涡流区o i 柏，故涡流区电场强度唯一。 电流源区，当盯l 两时电场强度e 唯一，0 ，= 0 时电场强度占的唯一性无法由此 式得到。 空气区电场强度的唯一性无法由此式得到。 3 2 泛函的推导和边界条件的证明 泛函：( 各向同性媒质) 齐次边界情况：。e ：d 或n x h ：d 或j 行。e 5 。”s 7 l n h = 0o ns 2 p 2 舌p 。v e v e dv + 号j 口p e d v + j c o i e j d v ( 3 - 1 ) f r e + 枷) 2 号p 。f v r e + 7 7 “2d v + 号1 e + 玎 1 2d v + jc o i ( e + 叩 j j d v 华北电力人学硕士学位论文 d f = l 。v e | v 咖+ j c o ol e - h d v + j 国i j d v 2 0 = d 一v er 肛h d s 十i h r v 。t 一甲e + jc o c y e + j c e d ) d v 2 0 由于体积分和面积分不同，上式被积函数应分别为0 。 即：v “一7 v e 十j o x r e + j c o s = o h 0 设面积分为k 。 当n e = 0 有h = 0 则k = 0 此边界须强加边界条件 当n h = 0 有h 0则k = d w 日h d s = o 此边界条件为自然边界条件 当i ln 。x e = 一口0 1 l s 1 k = f - h h d s + i 2 。v e h d s = 0 nh0o ns 2 k 2 l = 占， 则s 1 须强加边界条件，s 2 的边界条件为自然边界条件 非齐次边界情况：n e = a 或n 日= 6 或七n 。日e = ：6 a硎o n ；： 其中，丑jb 是位置的函数。 设m 为满足以上边界条件的任意函数，则令k = e m 满足齐次边界条件 则有v x 一v k + j c o a k = 一j 彬一v 。一v m 一，谢 代入3 - 1 式 p 2 与p 1 x f e m ) 习x ( e m ) d v + 七j e m ) ( e m ) d v + ，i r e m j r ，+ 面1v ，叫v m + c r m ) d v 展开并去掉与取极值无关的项。 f = 号i 一v e v e 咖+ 号- ，们i er e d v + j c oi e l ，咖 一j 卢一v e v m d v + 冒可一v m d v 应用矢量恒等式v a b = v r 爿bj + av b 得： f = 号阻。v e v x e d v + j f e e d v + j c o e ，j d v 一扣。v m 例e d s 当n x e = a 时面积分与取极值无关可消去，故x e = a 为强加边界条件： 当h = 西时f = 号阻一7 v e e d v + j c o c rf e e d v + j o ) f e - j d v + f b e a t s ， 故，l x h = b 为自然边界条件； 出j n n 。x e = ：a h o ns 1 f = 扎“勺x e 习xe ( 1 v + j l e edv+jh o ? ls 2 硝e 。d v 千婕e d s f 月= 6 。j 。 2 。 jj 毒3 故s 1 处须强加边界条件，s 2 处为自然边界条件。 综上：f 2 舌j 卢。v x 可e d v + - ；j c o c rj e e d v + j c o j 冒j d v + 曲出 ( 3 2 ) f e = 口o ns ， l 一h = b o ns ? 其中s 1 处须强加边界条件，s 2 处为自然边界条件。 不同媒质分界面： 巧f 中添加l h x 一v x e ，也d s l ，l x 弘一v e 2 h d s 上式2 r 月h ，- - n h ：) n d s 由于h 0 故月日，= 月h ，自动满足， 而，l e ，= n e ，由棱单元插值函数自动满足。 综上：媒质分界面条件由泛函自动满足。 华北电力大学硕士学位论文 3 3 泛函的离散 选用w h i t n e y 型长方体单元，单元模型及棱边编号如下：( 其中各棱边正方向同右侧直 角坐标系) 号 匕 其中n 7 = 【：f ，；，n i 露、e 7 = 陋；，e ；，e ；j 、 ；= 。，。：，。，。 、e ；= 陋。，e 。：，e 。，e 。 、n ；= i n n 心nr 3 ，n ，。】、 e j = b ，e 川e 川e ，。 、；= 【n ：，n ：，n 纠n ：。】、e ；= 魄，e 。：，e 。，e 。 、 n 1 = h 。( 一y ) ( p z ) 、 n 2 = h 。( y 一以) ( p z ) 、 n x 3 = 厅，( m y ) ( z 一叮) 、 n t = h 。( y 一”) ( z q ) 、n r l = h ，( p z ) ( 口一x ) 、n r2 = h ，( z 一目) ( 口一x ) 、n y 3 = h ，( p z ) ( x b ) 、 n r 42h ，( z g ) ( z b ) 、 n z l = h ：( a x ) ( m y ) 、 n z ：= h = ( z 一6 ) ( m y ) 、 z3 = h ：( 口一z ) ( y 一月) 、n z 4 = h ：( x 一6 ) ( y n ) 。 对3 - 2 式进行离散并对e 求导( 线性媒质) ，每个单元内有： 【p 。v v 7d v + j 们 r 7 咖归= 一j o ) _ f n j d v 一! 、r - b d s ( 3 - 3 ) 左侧第一积分项处理：令，7 = v x n 7 = n ii ，v n f _ ，v n ik j = | 警产等死等t 一警t 等z 一尝，1 、 f ，f t ：x n r x n t ： f “ f y x 疋。 f h o ， 巴， f e 屹 ：掣攀+ 等等、：碟：一譬掣、 o zo z 西钟o x ：譬警+ 譬譬、：吆：一譬攀、 。一百百+ i i “1 旷一百i 、 屹：等等+ 盟丝、：曙：一盟攀，o “ 砂砂xo x “” o z 却 具体积分结果参考文献 5 。 l o 华北电力大学硕士学位论文 3 4 电流源区处理 右侧第一积分项处理：( 仅对长方体型的电流源进行处理) 首先对电流源分区如下： 设电流为顺时针方向，设，= j x f + ，y ，+ j z k 、 j = | j l ，并且j = j x = j y ：j z 。 对于l 、3 、5 、7 区处理，以l 区的单元为例： - j 巩2 f 。，。d v 2 ，_ v ，d v = 导j 1 ，1 r 2 j t ， 丁= 等d ，1 ，1 ，1 7 斗 同理：3 区f n j d v = 一j t ，5 区f n j d v = 一旧，7 区j d v 对于2 、4 、6 、8 区处理，以2 区的单元( 如右图) 为例： 目先征副分盯保让此荚卑兀在x y 半回的截回为止万彤 公式的推导并且在剖分时这一点容易保证。 妒j d v 2 j i , n x d v j1 2 n y d v 其中j ，。咖= 西a 【l 2 1 ，2 r 、1 2 ，咖= 龛【1 ，1 2 ，2 】r 同理，4 区j n j d v 2 一jf ，d v 一，f ，n 。d v 其中f ，r 咖2 舍d ，1 1 2 ，2 】r 、c 。以咖= 舍 2 1 ，2 ，1 7 6 区j d v = l ，f ，n ，d v i ，f 。d v 其中f ，机咖= 盒盼l ，1 卜f 。以扯舍d ，1 2 ，2 7 8 区f - j d v = ，f ，。d v + jf 。n v d v 其中f ，以咖2 鲁n j 2 ，l ，2 r 、f 。y 咖= 舍 2 ，2 ，1 ，1 7 。 3 5 非线性处理( 牛顿拉夫逊法) | 8 i 2 卜 3 卜 s4 图33 电流源分区 这样做简化了 图3 - 4 电流源2 区 由于非线性磁阻率只出现在( 3 - 3 ) 式第一项积分中，即l 刃n 口e d v ，( 其中非线性 磁阻率表示为厂) 所以只处理此项。令厂= i 罗n v e d v 。 望：f v v e 里如+ p 旦v 可x e 咖 0 e o 0 e o a e 上式中右端第二项积分同线性情况，所以只处理第一项积分。 设磁通密度为b = b 。“b ，j + b 。k ，b = i 日。l + 1 b ，i + i b ：l = b 。否z + b ，百r + b 。百z ， 最2嚣嚣秀上瓦ob22b盟c3bkl瓦obxe beb2 bo e 商畿厩等e a x i 88 x 。ax i j 8 e x i。8 e x i 4 a x ! 华北电力人学硕士学位论文 = j 乏杀嚣( 否r 百o n x t 一百z 等 ，其中否x 表示b 。的共轭，在推导中应用了复变自 数的形式导数，参考 1 3 。 同理： 岳士1o y ( 1 也盟也盟) ，老io r ( i 否x 百0 n z i o e 2 c o bo bo zo xo e2 0 a bo b 巧r 丝o x ) ， ，； 。 lj， 7 la yj b ：，1 _ v n 7e 0 9 = ，去 ( 警夸尝毛 z + ( 警妒警。等妒警氐) 1 盟= 击h ，1 7 、一o n x = 击m - l1 卜百o n r = 西1 0 z 2 d o y 2 d2 d 1 ，l 】、 y 。 。 y 1 a z 7 。 1 尝= 击几- 卜百o n z = 瓦1 川7 、百0 n z = 击训。 y 与岳的求法参考 1 4 。一 一 ，ik ，k ，置，i 令彭= f v x n v x e 。0 。，咖。l 套：k k ：i t ，足k 。r zf k 。= j v 嘶m e 嚣7 咖 叫国户叱邶嘉7 咖 一，蠢嘉j ( 也百o n x b ：警) 蠢1 西 = 龛嚣( b 警一岛軎p 警五警 1 为了得到对称的系数矩阵，假设b 。b 。的相位相同。这一假设只影响牛顿一拉夫逊法的 迭代速度，并不影响解的精度“。引入假设上式可写成： j o r ( i 别警倒引o n x 、百0 n x 倒百o n xj 1 b 3 尘2 b 旦o b 。：警恻o 鳓n v 划百o n r 恻o 。n ：v 7 珏尘2 b 盟o b s o 一_ n z 旧l 等n 协o n z 制警 7 3 会剐吼i0 1 _ n x 峨i 警l 百o n r 也i 纠1 华北电力人学硕十学位论文 蝣聒= 尘a - - r ( i 玑 o 矿n zb i 尝n o i n x 恻o 万n x 翰魄气a b a b o v i ( i 引0 矿n zp 。i 百o n z ( 】b z ) l 百o n v 恻尝 1 另外的推导： k 。，v x n x i v x e 筹o y 咖 2 瓦1 瓦0 y 【b v 百0 n 。加z 警一警也等卜 2 五1 瓦0 y 7 m 百3 n x 百o n x r 刊| b z | ( 等等7 一等警7 卜l2 等等7 d y 。五1 面0 y1 酱蚶k 2 一鲁j b y l l 吲m o + 可d x d z | b z | 2 k 1 l 臣r r 2 r 2 v x n x i v x e u o l y y 7 d v 。去嚣m 。警也喾净警- - x 警卜 2 五1 瓦0 y b z l 百o n xi 0 n v 制刚百0 n x i 0 n v 倒2 百0 n x 0 苏n vr + | b z | | b 。i 百0 n xi 0 n y 7 d v 2 五1 瓦0 7 1 l 鲁吼忙朋t 一百d x f d v b y 忙抽t 2 一警阪陋3 + 手峨j 胛。1 k w 2 五1 而0 y7 d 。v 。d 。z 叫2 k 2 一了d y | b z | j b 。l n t 0 + 百d x d y 陬| 1 k 1 足。2 瓦1 瓦0 y1 百d z d x 附k 2 一了d z 刚叫n t 0 + 百d y d z 时k 1 k “瑙w r 。万1 瓦0 , vr | 鲁陬忙。f n t l 一鲁| b z | | 叫n t 2 - i d y 峨陬，+ 瓦d x d 叩z 。| | b z | 丁。 k y z 2 k z v , 一面1 瓦0 y l d 4 v d d 。z b z 陬f n t - 一鲁酬刚n t z 警刚2 k ，+ 鲁刚酬胛。 其中：k 1 、k 2 、k 3 的积分结果参考 1 。 j 20 0 2 111 1 一1 w t 。2 ：? 2 2：1w t - = i ! ，! 。- j 1 1j 华北电力人学硕士学位论文 n t 2 =n t 4 4 华北电力大学硕士学位论文 4 1 引言 第四章e 方法的程序实现 本文应用c + + 和v i s u a ll i s p 语言实现了上述有限元方法。其中c + + 用于主要计 算部分的实现，v i s u a ll i s p 用于前后处理部分的实现。二者的联系用数据文件实现。 在程序实现过程中对已有技术进行了改进和创新，提高了程序的总体效率。本章对 该部分内容进行介绍。 4 2 程序数据结构方面的改进 电磁场有限元方法通过对泛函方程或者加权余量方程的离散，归结为求解一组 实数或复数稀疏线性方程组。其中系数矩阵的稀疏性是有限元方程组的重要特征。 为了说明系数矩阵的稀疏程度，可定义矩阵的稀疏率为：s = 矩阵的非零元素数矩 阵的元素总数。又由有限元插值函数的局部性可知每个未知量和一定数目( 设为t ，t 与剖分所用单元相关，忽略边界情况) 的其他未知量相关。矩阵的稀疏率s 与方程 的阶数成反比，例如用一阶六面体棱单元对求解区域进行剖分，则t = 3 2 ，当阶数为 1 0 0 ，0 0 0 时s 仅为万分之三。若矩阵全部元素数的存储量为o ( n 2 ) ，若仅存储非零 元素，存储量仅为o ( n ) 。所以采用的稀疏矩阵存储技术对于减少系数矩阵的内存使 用总量非常重要。常用的稀疏矩阵存储技术大致可分为两类：数组存储技术( 如 1 j 中的行索引存储) 和链表存储技术( 如 2 、1 6 中的十字链表存储) 。两类技术各有优 缺点，但它们的优缺点具有互补性。所以若能综合应用两种存储技术则能够取长补 短，提高有限元算法的总体效率。以下将分为4 部分进行讨论，第4 2 节介绍上述 两类稀疏矩阵存储技术及其优缺点；第4 3 节介绍对十字链表的改进以及本文提出 的a c l 技术；第4 4 节介绍a c l 方法在以电场强度e 为变量的棱边有限元方法求解 时谐涡流场中的应用。第4 5 节为结论。本文算法是在w i n d o w s 操作系统以及c + + 环境下实现的。 4 2 1 常用稀疏矩阵存储技术简述 三维电磁场数值计算中，由一个一维或多维数组存储系数矩阵非零元素，再由 几个数组进行索引。不同的索引方式将对程序的执行效率产生很大影响。这类方法 的优点是分配单元连续，可以达到最小的内存占用量。缺点是固定的长度使得插入 和删除元素的操作耗费较多的机时，而这种两种操作对于有限元算法的系数矩阵形 成和边界处理过程非常重要。另外元素数组和索引数组的联系不如链表存储技术紧 密直观，增大了编程的难度。 华北电力大学硕十学位论文 链表存储技术晚于数组存储出现，但其灵活的特性使其在稀疏矩阵存储方面显 示出越来越大的生命力。链表存储一般采用内存动态分配技术，每个非零元素和它 的索引指针被定义成一个结构。添加和删除元素的操作都是指针操作，没有顺序的 限制，使得这两种操作非常自然灵活，可大大提高系数矩阵形成和边界处理过程的 效率。其缺点是链表占用的内存量要大于理论值，即用每个非零元素占用的字节数 ( 以下称为单元) 乘以非零元素总数并不等于非零元素占用的内存总量。因为内存动 态分配技术要占用一定的内存来存储分配信息，以及删除操作使得每次内存分配并 不连续使得单元之间形成无法再分配的内存空间。这样的内存空间约占非零元素总 内存占用量的2 5 。无形中增大了内存使用量且不易被察觉。而且，非零元素散布 在很大的内存空间中也会降低程序执行速度，因为现在的计算机普遍使用高速缓存 ( c a c h e ) ，c a c h e 的访问速度比主存储器快好几倍，但容量比主存储器小的多，c a c h e 只能按一定的规则更新存储主存储器的部分数据，当c p u 需要的数据在c a c h e 中( 称 为命中) ，读取速度比较快，若未命中，就要到主存储器中读数据，这时速度比较 慢。当数据的物理地址比较分散时，c a c h e 的命中率就会降低，必然影响程序的执 行速度。本文在第4 3 节将给出其定量分析，并提出a c l 技术弥补链表存储技术的 上述弊端。 4 ，2 2 改进十字链表技术 文献 2 中采用的十字链表法是一种双重循环表。稀疏矩阵中的每个非零元素 可用链表中的一个节点表示，这种节点由五个域组成，其中：行域( r o w ) 、列域( c 0 1 ) 、 值域( v a l ) 分别表示某非零元素所在的行号、列号和数值；向下域( d o w n ) 用以链接同 一列中表示下一个非零元素的节点，向右域( r i g h t ) 用以链接同一行中表示下个非 零元素的节点，如图4 1 示。于是稀疏矩阵中同一行的非零元素通过向右域链接成 一个带表头的循环链表，同一列的非零元素通过向下域链接成一个带表头的循环链 表( 行、列表头节点如图4 2 所示) ；另外文献 2 还对链表的操作( 包括插入、查找、 删除、求和、修改) 定义了统一的函数。这便形成了十字链表的基本结构。 对十字链表的改进： 在有限元过程中除了矩阵形成过程，其它过程只需卜- 号罴 j 产 孟h 行列表头结点之一，可在矩阵形成后释放掉行表头结点图4 - 1 非零元素节点 所占用的内存。另外对指定行( 列) 的查找操作对于矩阵r 合成和边界处理过程都非常重要，而表头结点应用链表丽f 上_ = 丽1 的形式，查找必须顺序进行，这样极大地降低了程序的 图4 2 行( 列) 表头节点 效率，使得上述两个过程的计算机时间可与方程求解的计算机时间相当。而且行f 列) 表 头节点的结构也过于复杂。所以要重新设计行( 列) 表头节点。在矩阵形成之前，矩阵的 阶数是已知的，这样就可动态分配两个指针数组作为行( 列) 表头数组。在矩阵合成之后 l6 华北电力人学硕士学位论文 可释放行表头数组。在进行边界处理时，由于行( 列) 消去操作使列表头数组产生了许多 空指针。这时可先再动态分配一个长度为当前矩阵阶数的指针数组指向矩阵，再释放掉 原列表头数组。这样在矩阵合成和边界处理过程中进行查找操作等于对数组进行索引。 相对链表结构的行( 列) 表头节点减少了大量的计算机时间，这使得矩阵合成和边界处理 过程的计算时间相对求解方程组的计算时间可忽略不计。改进后的十字链表结构如图 4 3 所示。 4 2 3 数组与改进十字链表混合的存储技术 在第4 2 1 节中已经叙述了链表存储技术的缺 即因为内存动态分配技术或删除操作使得每次内存 分配并不连续，这样就有可能在结构与结构之间形 成了无法再分配的内存单元。使得程序占用的内存 总量超过我们的计算值。具体数值可参见表1 。 a c l 技术可方便地解决上述问题。a c l 的主要 图4 - 3 改进的十字锛表结构 思想是对于程序中不需删除操作的单元( 称为永久单元) 先以数组的形式统一分配 内存，执行插入操作时，在数组的范围内进行分配；对于在边界处理过程中可将被 删除单元( 称为临