（地球探测与信息技术专业论文）代数多重网格法在直流电法有限元模拟中的应用.pdf

摘要 有限元数值模拟方法是解决地球物理复杂问题的有效途径之一。 有限元分析主要包括三个过程，即网格生成和优化、有限元离散代数 系统的形成以及离散代数系统的求解。其中，有限元离散代数系统的 求解是有限元分析中重要的一步，直接关系到数值计算结果的准确 性。因此寻求一种高效的求解方法，对直流电法有限元模拟具有重要 意义。 本文系统地分析了目前求解直流电法有限元模拟中的代数方程 组的常用迭代法，如c g 法、i c c g 法等。在对常用解法进行研究的 基础上，引入了具有高速收敛性的代数多重网格法( a m g ) ，对a m g 法的基本原理、收敛性和具体算法流程进行了详细的论证，并应用数 值算例对其计算效果进行了分析。通常情况下，直流电法限元模拟中 代数方程组通常具有大型稀疏特征，本文针对这一特性，研究了适用 于有限元代数方程组的a m g 法预处理共轭梯度法( a m g p c g ) 。 最后在m a t l a b 平台下编制了多种常用迭代法、代数多重网格法和代 数多重网格法预处理共轭梯度法的算法程序，并应用到实际模型计算 中，对各种算法的收敛速度和收敛效果进行了比较。 结果表明，a m g 法和a m g c g 法在求解直流电2 5 d 有限元模 拟中生成的方程组时，精度高、收敛速度快；求解效果优于其他迭代 法。 关键词代数多重网格( a m g ) ，代数多重网格预处理，直流电法 2 5 d 有限元模拟 a b s t r a c t t h ef i n i t ee l e m e n tn u m e r i c a ls i m u l a t i o nh a sb e c o m eo n eo ft h e e f f e c t i v e w a y st o s o l v e c o m p l e xg e o p h y s i c a lp r o b l e m s t h ef i n i t e e l e m e n tm e t h o dc o n s i s t so ft h r e em a i np r o c e s s e s ：m e s hg e n e r a t i o na n d o p t i m i z a t i o n 、f o r m a t i o no ft h ef i n i t ee l e m e n td i s c r e t es y s t e m a sw e l la s s o l v i n ga l g e b r a i ce q u a t i o n s a m o n gt h e m ，s o l v i n ga l g e b r a i ce q u a t i o n si s a ni m p o r t a n ts t e p ，w h i c hd i r e c t l yi n f l u e n c e st h ea c c u r a c yo ff i n a lr e s u l t s t h e r e f o r e ，i t ss i g n i f i c a n tt os e e ka ne f f i c i e n tm e t h o df o rs o l v i n gt h e e q u a t i o n so fg e o p h y s i c a lf i n i t ee l e m e n ts i m u l a t i o n t h i s p a p e ra n a l y s i s t h en o r m a lm e t h o d sf o r s o l v i n ga l g e b r a i c e q u a t i o n so ff i n i t ee l e m e n ts i m u l a t i o nd e t a i l e d l y ，s u c ha sc g 、i c c ge t c a f t e rr e s e a r c h i n gn o r m a li t e r a t i v em e t h o d s w ei n t r o d u c ear o b u s t l y c o n v e r g e n tm e t h o d ，w h i c hi sc a l l e da l g e b r a i cm u l t i g r i dm e t h o d ( a m g ) ， w ed e s c r i b et h eb a s i cp r i n c i p l e s 、a s t r i n g e n c ya n da l g o r i t h m so f a l g e b r a i c m u l t i g r i dm e t h o d ( a m g ) d e t a i l e d l y ，a n dw ea l s oa p p l i e da m gt o n u m e r i c a le x a m p l e st oa n a l y z et h ec o n v e r g e n c e u s u a l l y ，t h ea l g e b r a i c e q u a t i o n so fd i r e c tc u r r e n tf i n i t ee l e m e n ts i m u l a t i o ni sl a r g e s c a l ea n d s p a r s e ，a c c o r d i n gt ot h i sf e a t u r e ，w ep r o p o s ea na m gp r e c o n d i t i o n e dc g m e t h o dt h a ti ss u i t a b l ef o rf i n i t ee l e m e n tl i n e a re q u a t i o n f i n a l l y w e w r i t ec o r r e s p o n d i n gc o d e sf o rs e v e r a la l g o r i t h m sw i t hm a t l a b a n dw e c a l c u l a t es o m er e a lm o d a l so f2 5 dd i r e c tc u r r e n tf i n i t ee l e m e n ta n a l y s i s u s i n gt h e s ea l g o r i t h m st oc o m p a r et h e s em e t h o d s t h er e s u l t ss h o wt h a ta m gm e t h o da n da m g p c gm e t h o dh a v e v e r yf a s tc o n v e r g e n c es p e e da n dv e r y h i g ha c c u r a c yi ns o l v i n gt h e a l g e b r a i ce q u a t i o n so ff i n i t e e l e m e n ta n a l y s i s ，t h e s et w om e t h o d sa r e m o r ee f f e c t i v ec o m p a r e dt oo t h e ri t e r a t i v em e t h o d s k e y w o r d s ： a l g e b r a i cm u l t i g r i d m e t h o d ( a m g ) ， a m g p r e - c o n d i t i o n e r ，2 5 dd i r e c tc u r r e n tf i n i t ee l e m e n ts i m u l a t i o n 原创性声明 本人芦i 明，所呈交的学位论文足本人在导师指导下进行的形f 究 1 ：作及取得的研究成果。尽我所矢| l ，除j ，沦义f ，特别加以标注和敛谢 的地方外，论文中不包含其他人已经发表或撰! 写过的研究成果，也不 包含为获得中南人学或喜他单位的学位或证书而使用过的材料。与我 、- t - k 同工作的同志对奉研究所作的贡献均已在论文中作了明确的说明。 作者签名： 桦 r l , 廿j j ：五年月五日 学位论文版权使用授权书 本人丫解li 南大学有关保留、使用：学位论义的规定，即：学校 r 仃十义保留。学位沦义并楸据幽家或湖南省有火部i 、j 舰定送交学位论义， 允许学位论义被禽阅和f m 蒯；学校口j 。以公们j 学位沦义的全部或部分内 容，可以采j 十j 复印、缩印或其它丁 段保存学位论文。 司时授权，川习科 学技术信息研究所将小学位论义收录剑r f l 幽学位硷文伞文数据席， 并通过删络向科：会公众提供信，心、j j & 务。 作椭礼绎叶 蝴签名m i i , 仕j j ：埠f 上月霁i i 中南人学硕士学位论文 第一章绪论 第一章绪论 地球物理学中的场值计算方法总的归结起来有两类：即解析法和数值法。由 于电磁场理论的复杂性，在电法勘探中具有解析解的地电模型是比较少的，例如： 均匀半空问中的球体、柱体、或简单的层状介质等。稍加复杂的地电模型则没有 相应的解析解理论公式，因此数值计算就成为电磁法正演中的重要分支。其中最 常用的方法有两种：有限差分法和有限单元法。这两种方法相比，有限单元法更 适用于物性复杂分布的地球物理问题，而且其求解过程也比较规范化，所以有限 单元法在各个工程领域中都获得了越来越广泛的应用。 有限元离散代数系统的求解，是有限元分析中的重要环节，是影响有限元分 析整体效率的重要因素。有限元离散系统通常将一个偏微分方程( 组) 的求解转 化为一个大型稀疏线性代数系统的求解。随着研究问题深度的增加，构造出来的 线性方程组规模越来越大，且往往是病态的。由于方程组的病态性，微小的误差 也会导致结果的巨大偏差。而我们在实际测量和计算中又不可能避免误差，因此 研究如何更好的消除误差，取得较为精确的数值解，有很大的理论价值和应用价 值。由于实际问题的需要，迫切的需要发展高效率的计算方法，于是经过众多学 者的努力，产生了一种在多水平层次上消除误差的算法一代数多重网格法。 1 1 国内外研究现状 最初的有限元法是在5 0 年代弹性力学应用方面发展起来的。 o c z i e n k i e w i c z i l l 在力学，热学，流体力学等方面系统地给出了相应有限元方法 的理论分析，使得有限元方法形成了一个比较完整的体系。而有限单元法在地球 物理中的应用则是在7 0 年代才发展起来的。1 9 7 1 年c o g g o n 2 】首先将有限单元法 应用到地球物理上，基于其在力学等方面的应用推导出了适合地球物理领域的有 限元理论，导出了等价的变分方程并给出电磁场与激发极化异常模型的计算，最 终结果显示与h o h m a n n ( 1 9 7 0 ) 所用积分方程法所得到的结果一致，这为有限单 元法在地球物理学中的应用奠定了基础。然而，c o g g o n 所采用的计算网格缺乏 通用性，同时计算的速度与精度也达不到很高的水平。1 9 7 6 年蹦i 0 1 3 l 改进了点源 二维有限单元法，引入了一个通用性网格，使计算速度和精度都大为提高，从而 使有限单元法在电法j 下演中达到了一个适用的阶段。同年，r e d d y t 4 1 等人用有限 单元法实现了对三维电磁( e l e c t r o m a g n e t i c ，e m ) 问题的模拟。但是由于三维 e m 问题的复杂性( 不同于电阻率法中所求解的量为标量，e m 中为复杂的矢量) ， 中南人学硕+ 学位论文第一章绪论 所以在当时计算条件下，模拟的实现是以时间为代价的。为了提高计算的速度， 一种对于任意二维地质体( 源为线源，即2 d 问题) 的瞬变电磁响应新算法被提 出( g o l d m a n ，1 9 8 6 ) ，然而不足的是该算法不能应用于极化介质中。1 9 7 8 年， h m b b b y l 5 l 在c o g g o n 工作的基础上给出了轴对称体直流( d c ) 电阻率模型的有 限元计算，且该方法非常灵活，能够计算任何具有轴对称体的视电阻率，从而扩 展了有限元方法在电法勘探中的应用范围。有限元算法在电法中的应用的不断发 展也随着计算机的发展而不断完善起来，对于复杂的三维问题也有了相应的有限 元j 下演模拟结果( p r i d m o r e ，1 9 8 1 ) ，就计算时间上来讲要比积分方程更加经济， 而且对于电导率比值很大的非均匀地质体的模拟也比积分方程更加实用。 h o l c o m b 与j i r a c e k ( 1 9 8 4 ) 进行了稳定电流场的带地形的三维模型的视电阻率 计算。w a n n a m a k e r 与s t o d t ( 1 9 8 4 1 9 8 6 ) 、c h o u t e a l 与b o u c h a r d ( 1 9 8 8 ) 就大 地电磁法中地形起伏下异常体的电磁响应进行了二维有限元数值模拟，提出了基 于有限元数值解的大地电磁响应的起伏地形校正方法；g a n q u a nx i e ( 1 9 9 8 ) 及 其研究组提出了基于边界一区域的有限元法对大地电磁模型进行了相应的模拟； s t a l n a k e r 与e v e r e t t ( 2 0 0 1 ，2 0 0 2 ) 采用四面体线性单元对三维可控源电磁模 型问题进行了探索；2 0 0 5 年以来，d e m k o w i c z 与r a c h o w i c z 等人开展了稳定电流 场及电磁场的h p 型自适应有限元算法研究，对石油勘探中的直流电测井模型进 行了数值计算，得出了高精度的视电阻率结果。k e y ( 2 0 0 6 ) 开展了大地电磁二 维模型的自适应有限元计算，l i ( 2 0 0 6 ) 给出了海底可控源电磁法二维模型的初 步计算结果。 国内最初的有限元方法是由冯康提出，而其在电法勘探中的应用则是在7 0 年代引入并逐步发展起来的。最初，有限单元法多数应用于物探数字模拟或仅适 于对称体地电条件下视电阻率的计算。周熙襄，钟本善【6 l 等人将有限单元法应用 于任意二维地电断面和任意地形条件下视电阻率法的数值模拟计算中，并采用非 均匀网格剖分，最后结果表明对于绝大多数的地形、非均匀地质体，不同电阻率 分界面的地电断面均能较好的模拟，这为直流电法中有限元法的广泛应用奠定了 良好的基础。周熙襄【7 】等学者近一步将此方法应用于点源二维地质体( 2 5 d ) 的电 阻率法数值模拟中，并采用混合边界条件、三r 分解线性方程组等优化措施，使 2 5 d 电阻率法正演计算的速度和精度都比国外流行的l r i j o 方法和程序有了进 一步的提高。 在稳定电流场有限元模拟中，离散代数系统的求解是很关键的一步，最终会 影响视电阻率值得准确性。通常我们采用迭代法对方程组进行求解，较常用的方 法有c g 法、i c c g 法等。对于小规模问题以及病态程度不大的方程，这些方法通 常都能够准确求解。但是对于规模比较大、病态程度比较高的问题，上述方法的 2 中南人学硕士学位论文 第一章绪论 收敛速度往往很慢，其原因是在迭代过程中，误差的低频分量不能够被快速消除， 且随着迭代次数的增加，其收敛速度越来越慢，甚至于不能达到精确求解的目的， 最终影响有限元模拟的精度。因此寻找一种能够快速消除误差，加快收敛速度的 解法，对于稳定电流场的有限元模拟来说，具有重要意义。 代数多重网格法是在几何多重网格法的思想上建立起来的【8 】【9 】【i o l 。2 0 世纪 6 0 年代以来，大型稀疏线性系统的求解问题越来越多，对于实际计算中的一些 问题，传统的单层次的求解方法已经不能满足求解需要，经过大量的理论研究和 实践，人们提出了几何多重网格法，该方法在多级网格上消除误差，能够有效消 除迭代过程中误差的低频分量部分，达到快速收敛的目的。在几何网格法中，网 格是对应于实际的几何结构的，细网格由粗网格加密而成，迭代过程中的各算子 对应于特定的几何结构。因此应运几何多重网格法，必须要知道隐含的几何结构， 针对不同的问题需要编制不同的程序，对于很多问题，通常难于利用其复杂的几 何信息束定义多重网格法的各个分量，甚至难于找到合适的磨光迭代算子，这就 限制了多重网格法的推广使用。于是人们开始寻求一种单纯依赖于线性系统的多 重网格法。 代数多重网格法( a l g e b r a i cm u l t i g r i dm e t h o d ) 采用纯代数方法来划分粗细 网格，它不依赖于所求解问题的几何或者物理性质，仅仅是利用代数方程中的系 数矩阵性质来构造一套虚拟的网格，然后按照几何多重网格循环过程求解。因而 对于不同类型的应用问题，具有算法上通用性。 最初的利用系数矩阵自动生成虚拟网格的尝试，始于2 0 世纪8 0 年代，当时 将依赖于系数矩阵所生成的插值算子与g a l e k i n 法相结合，以提高算法的收敛速 度。第一个有关于a m g 法的程序是在2 0 世纪8 0 年代出现的，此后很多年罩， a m g 方法几乎没有新的进展，直至2 0 世纪9 0 年代中期，随着几何多重网格法 越来越不适用于复杂性问题，其应用的局限性越来越明显，人们丌始把目光投向 代数多重网格方法。经过k s t u b e n 和j r u g e 8 】【9 】等人的努力，形成了目前被称之 为“标准a m g ”算法的代数多重网格算法。标准a m g 方法的粗化策略是通过设 定强弱联通的阈值，对系数矩阵的非零元素进行判断，将其对应的虚拟节点划分 为粗网格节点或细网格节点，由此实现纯代数水平的粗细网格划分，从而在此基 础上构造多重网格算法所需要的算子。其中粗细网格的划分过程通常要占用大量 的时间，对于大型稀疏矩阵较为实用，对具有较多非零元素的矩阵，在时问丌销 和存储量方面的要求很高。 2 0 世纪八九十年代，代数多重网格法在“标准”算法的基础上得到了进一 步的发展。在此期间，出现了基于光滑聚类法的a m g 法( a m gb a s e do ns m o o t h e d a g g r e g a t i o n ) l ，以及基于部分几何信息的a m g 法( a m g e ) 1 1 2 】。 中南大学硕十学位论文 第一章绪论 此外还有很多有关于求解某种特定类型方程组的a m g 方法【3 l 【1 4 l 【1 5 】【1 6 1 ，这 些方法的求解过程是一致的，区别在于粗细网格划分方法的不同。 “高效性”和“快速收敛性是a m g 方法始终追求的目标，近年来a m g 方 法主要朝着两个方向发展，一是发展新的适合于不同问题的粗化、插值方法；二 是进行并行计算以及分布式计算0 7 1 1 8 l 【1 9 】【2 0 l ，开发并行a m g 以及分布式a m g 计算方法，提高a m g 方法的求解效率，比较有代表性的并行算法有c l j p 2 ， b o o m e ra m g 2 2 1 。 目前已有的各种针对不同问题的a m g 方法，有着不同的运行效率和数值稳 定性，但是没有一种适用于任意问题的a m g 解法存在。在很多情况下，其求解 效率往往会大打折扣。 为了增强多重网格法的收敛性，9 0 年代以来人们越来越多的将多重网格法 作为c g 法的预处理方法【8 1 ，而不是单纯的使用多重网格法进行求解。大量的实 验证明，利用a m g 方法构造c g 法的预处理矩阵，能够达到快速收敛的效果。 有关于a m g 做预处理的c g 法，我们将在后续章节详细介绍。 1 2 本论文的主要内容 本文第二章介绍了稳定电流场的电位方程、边值问题，以及稳定电流场2 5 d 有限元分析的基本理论和实现步骤。 第三章介绍了目前求解大型方程组的常用迭代解法，包括y a c o b i 迭代， g u a s s s e i d e r 迭代法等，并对病念问题以及预处理方法的理论进行性了描述，详 细推导了预处理共轭梯度法的计算过程。 在第四章中，我们介绍了多重网格法的基本思想，给出了代数多重网格法的 求解过程，分析了多重网格法的收敛性，并详细介绍了“标准a m g ”网格粗化算 法以及基于光滑聚类法的网格粗化算法。给出两个数值算例，对多重网格法的求 解效果进行分析。针对直流电有限元分析中的大型稀疏方程组求解，引入了代数 多重网格法做预处理的共轭梯度算法。 在第五章中，我们运用代数多重网格法和代数多重网格做预处理的c g 法求 解了直流电有限元分析中的一些典型模型，并将代数多重网格法和代数多重网格 做预处理的c g 法的求解效果同常用迭代法的求解效果进行了比较。 4 中南大学硕士学位论文 第二章稳定电流场2 5 d 有限元模拟 第二章稳定电流场2 5 d 有限元模拟 2 12 5 维稳定电流场的边值问题 2 1 1 稳定电流场边值问题 在线性各向同性介质中，描述电磁场的m a x w e l l 方程组为乜引： ( v 日) = j + 百0 1 ) ( 2 - 1 ) v 西：一竽 ( 2 2 ) 西 、。 v ( d ) = q( 2 3 ) v ( t h ) = 0( 2 4 ) 其中，日为磁场强度，_ ，为传导电流强度，d 为电位移矢量，e 为电场强 度，b 为磁感应强度，q 为体电荷密度。 在电磁场中还应考虑介质与电场和磁场的关系。其关系方程式为： d = c e ，b = t h ，j = o e 其中o r 为电导率，为磁导率，s 为介电常数。 对于稳念电流场的情况可知： v e = 0 ( 2 5 ) d a ( 2 5 ) 式可知，存在一个标量势函数”使得 占= 一v u ( 2 6 ) 对欧姆定律j = e r e 两边取散度并将式( 2 6 ) 代入可得： v ( o v u ) + v ，= 0 ( 2 7 ) 10 萑q 利用万函数性质： 烈锄d q2 t 圭 q 如图2 1 r q 分别为模型区域及 可得：对于含源区域( 如为单一的点源i ) 有心4 1 ： 中南大学硕士学位论文 第二章稳定电流场2 5 d 有限元模拟 其中a = a ( x o ，y o ，z 0 ，) 为源的位置，将上式代入式( 2 - 7 ) 中可得： v ( d v 甜) = - 2 1 8 ( a ) ( 2 8 ) 式( 2 - 8 ) 式即为电位的控制微分方程。 在直流电法( d i r e c tc u r r e n tr e s i s t i v i t ym e t h o d ) 中，加上边界条件可得 稳定电流场中的边值问题2 叭舶1 。 图2 - 1 三维模型 如图2 一l ，假定地面水平，坐标系建立如上所示。其中r 。，r ，r ，分别为截断 边界、两介质交界面、地一空边界；碍，为两介质交界面处的外法向方向矢量。 控制微分方程： v ( o r 甜) ：导仃搴+ i 0 盯i o u + i 0 盯j o u ：- 2 2 6 ( a ) q ( 2 9 ) 7 魂苏却却瑟瑟 、 地一空交界面： 譬l = 0 e f 。 ( 2 - 1 0 ) 两介质交界面： 塑=一0201 堕1 7 01 = 一1 ，no n o n 。,o n 2 截断边界上混合边界条件： o _ u + c o s ( r , n ) u ：0 r 。 o nr 6 ( 2 - 1 1 ) ( 2 - 1 2 ) 中南人学硕士学位论文 第二章稳定电流场2 5 d 有限元模拟 其中，为点源到边界的距离：，= ( x 一硝) 2 + ( y n ) 2 + ( z - z 一) 2 ，n 为截断边界 的外法向方向矢量。( 2 1 2 ) 式推导如下： 在所取区域足够大时，无穷边界上的电位值应该可以看作是一次场 u = f 。( 2 - 1 3 ) 对上式进行截断边界的外法向方向上的求导： 学：掣 ( 2 - 1 4 ) 锄锄 、。 转化为对，的求导： c d 一 詈= 舌一吾c o s 玎)a na rr 2 一 对于双点电源( ，万( 爿) ，一i s ( b ) ) 的情况，( 2 1 2 ) 式变成： ( 2 - 1 5 ) 冬+ 忐( 尘c o s ( _ ，刀) 一生c o s ( ，胛) ) = o ( 2 1 6 ) 咖【一r aj 屹 其中，_ ，r a 为点源到求解点的距离，对于多点源情况可以此类推。 2 1 2 稳定电流场2 5 维边值问题 假设地质体电性分布沿z 方向不变，在单点电源供电时，则位函数u 满足的 方程为幽3 3 ： 昙( 仃鱼唾趔) + 昙p 塑唾幽) + 昙p 塑唾塑) ：峨, y a , z a ) ( 2 - 1 7 ) 反、良 7 咖、西 7 岔、岔 各个边界条件与三维的边界条件( 2 一1 0 ) 一( 2 1 2 ) 相同。 对方程( 2 一1 7 ) 两边对z 进行傅氏变换。假设位函数对称于x o y 平面，则傅氏 变换可以化为余弦变换，积分区间简化到( 0 ，0 0 ) ，从而有 u ( x ，y ，k ) = “( x ，y ，z ) c o s ( k z ) d z ( 2 - 1 8 ) 记f ( u 户u ，对( 2 - 1 7 ) 式进行傅氏变换： 由于兰仃掣，晏仃掣分别是对x ，y 求偏微分与z 无关，故可以将 u 再 6 0 f , t o u y 积分与偏微分符号顺序互换，其傅氏变换为 7 中南人学硕十学位论文第二章稳定电流场2 5 d 有限元模拟 f ( 旦盯刿) = 旦盯o u ( x , y , k ) ( 2 - 1 9 ) 、融苏 7 苏 缸 、 f ( 晏仃掣掣) ：晏盯_ o u ( x , y , k ) ( 2 - 2 0 ) o ) o ) 哪 对于晏伊o u ( x ，, y , z ) y 页，由于是对z 求偏微分故不能将积分与偏微分顺序进行调 换。利用分部积分及条件：譬i：o ，i o u i ：o ，u ( o o ) ：o 可得 c ! z l ：；o出i ：。 ，( 兰盯掣) ：一七2 0 - u ( w ，七) ( 2 - 2 1 ) 由万函数的性质： ( 矽( 儿) 万( 乃) c 。s：昙万(_)万(n)(2-22)(kz)dzi 万( 矽( 儿) 万( 乃) c o s= 去万( _ ) 万( n ： z 故( 2 - 1 7 ) 式经傅氏变换后为： v ( 田u ) 一k 2 c r u = - 1 8 ( a )( 2 - 2 3 ) 其中，u = u ( x ，y ，k ) 。这个方程中只含有x ，y 坐标，所以为二维方程。方程中k 为参数，关于参数k 的选择，见文献 2 5 。 对边界条件( 2 一1 0 ) 进行傅氏变换：由于偏微分是针对法向方向，与z 无关， 故可以交换积分与微分的顺序。经傅氏变换后如下： 掣i ：0 ( 2 - 0 2 4 ) i l 2 u l lh 对边界条件( 2 1 1 ) 进行傅氏变换与上述同理可得： q 警= - - o 2 警f ，( 2 - 2 5 )q ? i ， 对边界条件( 2 1 2 ) 进行傅氏变换，式左边第一项与上同： ，( 盯宇) ：盯掣 ( 2 - 2 6 ) 、锄7锄 、7 由公式“：上可得： r ( x ，y ，z ) 肚f ( 南) = j c o 南c o s ( 舭= c k o ( k r ) ( 2 - 2 7 ) 其中r 是通过坐标原点且与走向垂直的截面内f 。上的点到电源点的距离。据 8 中南人学硕士学位论文第二章稳定电流场2 5 d 有限元模拟 旱k ( x ) = 一墨( x ) ( 2 2 8 ) 则( 2 - 1 2 ) 式第二项变为： 型+ 后梨c 0 妣刀渺：o r 。 ( 2 2 9 ) o n k ( 扫) 一7 。 、 7 同理，式( 2 一1 3 ) 经傅氏变换后为： 盯筹+ 盯币丽k u k ( 帆) c 。s ( ) 一k ( ) c o s ( 训= 。( 2 铷) 其中厶，分别是通过坐标原点且与走向垂直的截面内观测点至电源点a 、b 的 距离。对于多点源情况依此类推。 2 2 2 5 维边值问题的有限元方法 2 2 1 基于伽辽金加权余值法的有限元方程 单点源情况下，2 5 维边值问题的伽辽金加权余值表达式为他引： j = ，i s ( a ) 6 u d q 一，k2 a u s u d q 一v ( s u ) v ( a u ) d n nnn r ( 盯! 王墨掣u ) 万u d r j = 0 对于双电极供电情况( i s ( a ) ，一舾( b ) ) ，上式变为： j = f i ( 6 ( a ) - 6 ( b ) ) 6 u d n f v ( s u ) v ( 仃u 矽q p 2 c r u d u d q qnn ( 2 3 1 ) 【k ( k r , ) c o s ( r _ ，”) 一k 。( k r , ) c o s ( r 8 n ) s u d r ( 2 - 3 2 ) 9 一急一帆 一k 盯 0 fj = 中南大学硕士学位论文第二章稳定电流场2 5 d 有限元模拟 图2 - 2区域q 离散为各个单元e q 一 将区域q 剖分为有限个小单元e ( 如图2 2 ) 。在单元e 上，对( 2 2 5 ) 式进行 积分，然后对各单元求和。单元上的电位可由节点上的电位插值得到： u = m 吖= n t 玑 ( 2 3 3 ) 对上式两端求变分： 8 u = n 1 6 u , ( 2 - 3 4 ) 对式( 2 - 3 4 ) 两端求梯度： v ( 万u ) = 刚1 万玑 ( 2 - 3 5 ) 代入( 2 3 1 ) 中可得 i i s ( a 矽t r 0 1n d q p u ? v n o v n ? u c d q l k 2 a 6 u ? n n l u o a n 二r ! 万玑tc 盯竺量曼乎t 玑，d r ：。2 ( 2 - 3 6 ) 将上式写作 万玑t k 。玑+ 万阢t k 2 。玑+ 万玑t k 3 。玑= 万玑t 只 ( 2 3 7 ) 其中 k i 。= i v n o v nd q ( 2 - 3 8 ) 五。= f k 2 0 - n n t d f ( 2 3 9 ) k ，p 等删 ( 2 - 4 。) 只= ，h 罗( a ) n d t a ( 2 - 4 1 ) i o 中南大学硕士学位论文第二章稳定电流场2 5 d 有限元模拟 将( 2 - 3 7 ) 扩展成有全体节点组成的矩阵或列阵，然后将各单元相加，得 8 u t k u = 万u t p( 2 - 4 2 ) 从( 2 - 4 2 ) 式中消去8 u t 得到基于伽辽金法的有限元法方程： ku=p(2-43) 对于双点电源的情况推导相同故在此只给出结果： 万仉t k i 。玑+ 6 阢t k 2 。玑+ 万玑t 墨。阢= s u e t 只 ( 2 4 4 ) 其中： 五。= v n l d q ( 2 4 5 ) k 2 。= ，k 2 a n n t d f ( 2 - 4 6 ) f ：e 墨e ，! 盯瓦丽万k j 而k ( 饥) c 。s ( 一k ( ) c 。s ( 刀) m 7 d r ( 2 - 4 7 ) 只= j ( ，万( 彳) 一1 6 ( b ) ) n d n ( 2 - 4 8 ) 同单一点源相同，将式( 2 4 4 ) 扩展成有全体节点组成的矩阵或列阵，然后将 各单元相加，得 8 u - k u = 6 u 1p ( 2 - 4 9 ) 消去8 u t 得 置u = p ( 2 - 5 0 ) 对于多点电源情况同理。 2 2 22 5 维有限元单元分析 1 ) 单元形函数2 8 1 冽 本文中采用的是三角形网格剖分，采用线性插值。如图2 3 对于任意一单元 e ( i ，j ，1 1 1 ) ，电位可以表示为： 甜= 甜+ 砂+ c ( 2 5 1 ) 中南大学硕士学位论文 第二章稳定电流场2 5 d 有限元模拟 y b j ，y j 、) 图2 - 3 任意单元在总体坐标系下的表示 将三个节点的电位及坐标值代入可得 j = a x j + b y i + c 材，= a x , + b y , - i - c u 。= a x m + 饥+ c 由上三式可得 口= 去【( 乃一虼) + ( 一片) “，+ ( 片一乃) 】- - 去a ( a a u i 4 r a j u j 4 - 甜。) 6 = 云1 【( 一誓，) 甜，+ ( 一) 甜，+ ( 一f ，) 甜，】- c b , 甜，+ 乞u j + 6 肿“。) c = 去 ( 一虼一只) + ( 片一虼) 坼+ ( 乃一一片) 】= 去( c ，+ q 叶+ ) 其中 n | 2y i y m ，b t = x m x j c i = x y m x m v | a | = y m y i ，b j = x i x m ，c j = x j , 一x l y m a m = y i yj ，b m = x j x j ，c m = x , y l x j y i = 三( 口屯一a j b ，) 将( 2 5 5 ) ( 2 6 1 ) 代入( 2 5 2 ) 中得： ( 2 - 5 2 ) ( 2 5 3 ) ( 2 - 5 4 ) ( 2 5 5 ) ( 2 - 5 6 ) ( 2 - 5 7 ) ( 2 - 5 8 ) ( 2 5 9 ) ( 2 6 0 ) ( 2 6 1 ) 甜= 去【( q x + 6 ，少+ q ) + ( q 工+ 乞y + 巳) 材，+ ( x + k y + ) 】 ( 2 6 2 ) = n i + nj u j + n 一m = l p i + l p j + l m 中南大学硕士学位论文第二章稳定电流场2 5 d 有限元模拟 其中，m ，m 称为单元p 似圳的形函数。厶，l ，l 为面积坐标。形函数具有 如下的性质： j = 铝客乡 单元节点的形函数性质如图2 - 4 ： 图2 - 4 单元形函数的性质 ( 2 6 3 ) 2 ) 单元矩阵分析 对( 2 - 3 7 ) 一( 2 - 3 9 ) 进行单元内的积分则 镌= 少警警+ 警警，= 盯- - 去( a p a q + b p b q 胁州。m ，c 2 删 k h ：k 2 盯 n i n i n l n jn | n m n ，n n jn ，n ，i = a k 2 盯 n m n n m n | n m n m 1ll 6 1 21 2 lll 1 2 61 2 l l1 1 2 1 26 ( 2 - 6 5 ) 对于式( 2 - 4 1 ) 来说，假设边i j 落在无穷边界上如图2 - 5 ( 其中q 为求解场区，r 为 求解区的边界) 。 图2 - 5 节点位于边界上情况示意图 1 3 中南人学硕十学位论文第二章稳定电流场2 5 d 有限元模拟 落在截断边界上，节点离电源距离足诋故可将等c o s 循作常数。 此时r 可取点电源到i j 边的中点的距离。因为为边界积分则此时的积分项中m ， m ，虬不再为原来所足义的形函数，为此，在i j 上应用等参单兀进行线性插 值。其中 m ：监2 ( 2 - 6 6 ) n j = 了1 + 4 ( 2 - 6 7 ) i | x = m + m ( 2 - 6 8 ) y = n i y i + ni y j d r = 忑夏丽= 丢d 善 其中，为i j 边的边长。 式( 2 - 4 0 ) 积分为 ( 2 6 9 ) ( 2 - 7 0 ) _ ic r k j 鼍k 铲( k rm 七型2 k 掣f l i i 麓籀 川， r 。o )o ( i 呻 l ，i 。 、。 扩展后 孟h ：盯七竺! 坐! ! 竺! ! ! ：叫l 6 i l 3 k 。( k r ) o 据万函数的性质，式( 2 - 4 0 ) 的积分可表示为： ，万( 彳) n d q = 乏1 ，( 彳) q 一 ( 2 - 7 2 ) ( 2 7 3 ) 即只有点电源处才有贡献。 3 ) 总体刚度矩阵的合成 对( 2 3 7 ) 中的各项进行扩展，并将各单元的贡献进行叠加可得最后形成总体 刚度矩阵。 假设最后剖分的节点总数为n ，如对任意三角形e ( i ，j ，m ) 1 4 中南人学硕士学位论文第二章稳定电流场2 5 d 有限元模拟 u 。= 【v 】r = 【o 弘】7 k 、k 毛e2 悟c 廿 l k k 。k k ，k p 乃2 。2 p 屯e2 k r z r 铲 h 扫。2 。2 。 卜铲k 3 。 露，e2 k c 矿矿 【- k ，k 3 8 k 。3 。 p = 【o 昙，0 】 则式( 2 2 7 ) 中有限元法方程矩阵可表示为： 眉：y 后，+ 七j e + 后， _ 一 j _ l 其中l i e 为模型剖分的单元总数。 1 5 ( 2 7 4 ) ( 2 - 7 5 ) ( 2 7 6 ) ( 2 7 7 ) ( 2 - 7 8 ) ( 2 7 9 ) 柳 研 埘 1。l 1。上1。o。l上 中南人学硕十学位论文 第三章求解有限元方程组的常用方法 第二章求解有限元方程组的常用方法 我们通常采用迭代法对方程( 2 - 4 3 ) 进行求解。经典的迭代解法有【3 0 】【3 1 1 1 3 2 】： j a c o b i 迭代法、g u a s s s e i d e l 迭代法、超松弛迭代法( s o r ) 、共轭梯度法( c g ) 等，本章将详细介绍以上经典迭代方法。 3 1 迭代法概述 迭代法是通过逐步迭代逼近精确解的，具有程序设计简单，适于自动计算， 是求解工程技术中产生的大型稀疏方程组的有效工具。 迭代法的一般形式： 对于方程血= b ( 3 一1 ) 将4 表示为矩阵m 与的差： a = m - n ( 3 - 2 ) 其中m 为可选择的非奇异矩阵，且使m x = d 容易求解，一般选择为a 的某种 近似，称m 为分裂矩阵。 于是求解线性方程组( 3 1 ) 化为求解m x = n x + b ，即求解 沪m 1 n x + m 1 b 若记 曰爿玎1 ，户竹1 6 则( 3 1 ) 可写为如下等价形式； x = b x + f ( 3 3 ) 由此构造向量序列 讧味。+ 1 钮x 。可 ( k - o ，l ，) 若缸。) 收敛至某个向量x 宰，则可得向量x 幸就是所求方程组a x = b 的准确 解。 迭代法收敛的充分必要条件是p ( b ) o 为可选择的松弛因子。 缈 s o r 迭代法的迭代矩阵为啪1 ： 瓯= 帆一心= ( d - c o l ) 。1 ( u + ( 1 - w ) d )( 3 - 9 ) 迭代格式为： “) = ( d c o l ) 一1 ( ( 1 一a , d ) + o j u ) x 0 + e o ( d 一彩工) 一bk = l ，2 ，( 3 1 0 ) 当缈= l 时，s o r 迭代法就是g u a s s s e i d e l 迭代法，因此，适当的选择松弛因 子f _ o 可使得s o r 迭代法比g u a s s s e i d e l 迭代法具有更快的收敛速度。 只有当松弛因子0 a j 2 时，s o r 迭代法收敛。 3 3 共轭梯度法 五十年代初，h e s t e n e s 和s t i e f e l 【19 5 2 提出了一种解线性代数方程组的新迭 代法，即为“共轭梯度法”( c g 法) 。实际上，c g 法是建立在最速下降法基础上 的，可看作是最速下降法的一种修正。下面我们先介绍最速下降法，并导出共轭 梯度算法。 3 3 1 最速下降法 在讨论最速下降法之前，首先讨论下降算法的一般思想【3 3 1 。考虑到某一迭 代点6 ( ，要产生新的点扩) ，使 m ( b l k + o ) o )( 3 11 ) 式中，p 为一向量，称为方向，f i 为一整数，称为步长。显然要使( 6 ( ) 0 ，置滓o ； ( 2 ) 计算搜索方向a = = 6 一如。 ( 3 ) 确定搜索步长f = 揣。 ( 4 ) 令焉+ ，= + f ，尼，若忙+ 一圳： 0 ，令k - o ； ( 2 ) 如果恢0 占，停止计算，输出工为近似解，终止算法，否则下行； ( 3 ) 计算= ( 4 ) 计算，( 川) = 6 一血( 川) ，屏= ( 5 ) 置k = k + 1 ，转( 2 ) 。 酽荔可 x 似+ 1 ) = x 似j + ，肛) p ( 。) ； ，( 川) = r ( k + l + 屏p ( 。) 理论上讲，对于n 阶线性方程组，利用共轭梯度法只需通过n 次迭代就可以 达到原方程的精确解，但实际上，由于计算机舍入误差的存在和矩阵4 的一些 病念特性，利用共轭梯度法迭代n 次后，通常得不到方程的精确解，以此，常将 共轭梯度法与预处理方法相结合，以改善方程性念，加快收敛速度。 3 4 预处理共轭梯度法 本节介绍方程的病态性质，以及改善方程性态性的常用方法，给出预处理共 轭梯度算法。 3 4 1 病态问题与条件数 在数值实验中，我们发现，对于某些方程组