（地球探测与信息技术专业论文）带地形直流电测深二维正演.pdf

f o r w a r do ft w od i m e n s i o nd i r c e tc u r r e n tr e s i s t i v i t y t o r w a r t lo lt w o - 1 1 i m e ni r c e tn r r e n tl k e s l s t l v l t y s o u n d i n g w i t ht o p o g r a p h y m a jo r ：g e o p h y s i c a lp r o s p e c t i n ga n d i n f o r m a t i o nt e c h n o l o g y d i r e c t i o no fs t u d y ：n u m e r i c a lm o d e l i n ga n di n v e r s ei m a g i n g g r a d u a t es t u d e n t ：l i a n g f a n g m i n s u p e r v i s o r ：p r o f o u d o n g x i n c o l l e g eo f e a r t hs c i e n c e s g u i l i nu n i v e r s i t yo ft e c h n o l o g y s e p t e m b e r , 2 0 0 8t oa p r i l ，2 0 10 8m 0叭7m7，m呲mmim 7 1y 研究生学位论文独创性声明和版权使用授权书 独创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作及取得的研究成 果。据我所知，除了文中特另t l d n 以标注和致谢的地方外，论文中不包含他人已经发 表或撰写过的研究成果，也不包含为获得其它教育机构的学位或证书而使用过的材 料。对论文的完成提供过帮助的有关人员已在论文中作了明确的说明并表示谢意。 学位论文作者( 签字) ：盈茑刍竺 签字同期：2 翌f 殳：垒：z 理 学位论文版权使用授权书 本学位论文作者完全了解( 学校) 有关保留、使用学位论文的规定，有权保留并 向国家有关部门或机构送交论文的印刷本和电子版本，允许论文被查阅和借阅。本 人授权( 学校) 可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采 用影印、缩印或扫描等复制手段保存、汇编学位论文。同时授权中国科学技术信息 研究所将本学位论文收录到中国学位论文全文数据库，并通过网络向社会公众 提供信息服务。( 保密的学位论文在解密后适用本授权书) 本论文是否保密：是 不 如需保密，保密期限为： 学位论文作者签名：凛荔锄 签字日期：乙pf o 年月加日 乡日 卜细 盆月 乏 彩砗 ：刎 字乳 签 刚岍 f， c 一 导 签 摘要 正演是电法反演问题和观测资料解释的重要基础，地球物理正演的数值算法，种类很多，有 限单元法由丁适用丁物性复杂分布的地球物理问题，且解题过程具有规范性从而成为正演模拟的 主要方法之一。 本文应用有限单元法进行点源_ 二维电场止演模拟计算，在j e 演过程中采川三角单元剖分网格， 有利丁对复杂地形的模拟，从而在往后的实测数据在反演时不需要进行地形改止；电导率( 口) 和 电何在相邻节点间线性变化，使得电导率和电位与实际地质构造更接近、与实际地电场变化规律 更易吻合；采用双网格系统，用人网格反映地下电性变化和测深点的位置，小网格进行实际的有 限元数值模拟计算；在5 个波数的有限元点源二维直流电测深总电位法模拟的基础上，提出一种 简单的对电位进行修止的方法。利j f j 均匀人地的理论值和有限元模拟值的比值作为修止系数，以 此来消除总电位法中电源点的影响，并且提高了在有限个波数的情况卜人极距模拟的准确性。在 点源二维电场止演的基础上，将一定地形条件下的二维地电模型的视电阻率响应，分解为畸变场 和止常场二部分进行分析，采川比值法消除由地形引起的畸变场，从而获得水平地形条什卜的视 电阻率值。 并通过几个模型计算，证明正演模拟方法的正确性以及地形校正的效果。 关键词：正演；有限单元法；数值模拟：直流电测深：地形校正 a b s t r a c t ： f o r w a r di st h eb a s eo fi n v e r s i o na n dt h ee x p l a n a t i o no fo b s e r v e dd a t a t h e r ea r em a n y k i n d so fn u m e r i c a lu s e df o rg e o p h y s i c a lf o r w a r d p r o b l e m w i t hn o r m a lt u r n a r o u n d p r o c e d u r ea n ds u i t e df o rt h ep r o b l e mw i t hc o m p l e xe l e c t r im e d i ad i s t r i b u t i o n ，f i n i t ee l e m e n t m e t h o d b e c o m e st h eo n eo f p r i n c i p l em e t h o d so ff o r w a r d i nt h i sp a p e r , a c c o r dt ot h ef i n i t ee l e m e n t sm e t h o df o r w a r dm o d e l i n go fp o i n t s o u r c e e l e c t r i c a lf i e l do nat w o d i m e n s i o ns t r u c t u r e ，a p p l y i n gf e mf o rf o r w a r ds i m u l a t i o n i n f o r w a r do p e r a t i o n ，t h et r i a n g u l a rc e l li se m p l o y e df o rt h ec o n v e r s i o nc a l c u l a t i o no ft h e2 一d f e m ，s ot h a tn ot o p o g r a p h i ca d j u s t m e n ti sm a d eb e f o r et h ei n v e r s i o no ft h em e a s u r e dd a t a a n dw i t he l e c t r i c a lc o n d u c t i v i t y ( g ) a n dt h ee l e c t r i cp o t e n t i a li nt h et r i a n g u l a rw i t hl i n e a r c h a n g e ，w h i c ha l l o w s e l e c t r i c a l c o n d u c t i v i t ya n dp o t e n t i a lg e tc l o s e r t ot h ea c t u a l g e o l o g i c a ls t r u c t u r e ，a n dt h ea c t u a lc h a n g ei ne l e c t r i cf i e l di sm o r ea c c o r d a n t d u a l g r i di s u s e df o rt h es i m u l a t i o nt h a tt h ec o a r s e 鲥dr e f l e c te l e c t r i c a l c o n d u c t i v i t yd i s t r i b u t i n g u n d e r g r o u n da n dt h es o u n d i n gp o i n t s ，w h i l et h ed e n s eg r i du s e df o rf o r w a r dc a l c u l a t i o n b a s e do nt h et w o - d i m e n s i o nt o t a lp o t e n t i a lf e m ( f i n i t ee l e m e n tm e t h o d ) s i m u l a t i o nu s i n g 5w a v en u m b e r so fp o i n ts o u r c ed i r e c tc u r r e n tr e s i s t i v i t ys o u n d i n g ，as i m p l em e t h o dt o c o r r e c tt h ep o t e n t i a lf r o mf e mi sp u tf o r w a r d u s i n gt h er a t i ob e t w e e nt h e o r e t i cu n i f o r m e a r t hp o t e n t i a la n df e mv a l u e sa sc o r r e c t i o nc o e 伍c i e n t st oe l i m i n a t et h ei m p a c to fp o w e r p i o n ti nt h et o t a lp o t e n t i a lm e t h o da n di m p r o v et h ep r e c i s i o no ff e mv a l u e si nt h ec a s eo f l i m i t e dw a v en u m b e r sa n dl a r g ee l e c t r o d e sd i s t a n c e b a s e do nt h ef i n i t ee l e m e n t sm e t h o d f o r w a r dm o d e l i n g ，w ed i v i d e dt h ea p p a r e n tr e s i s t i v i t yr e s p o n s eo f2 dm o d e li n t ot w op a r t s ： n o r m a lp a r ta n da b n o r m a lp a r t ，a n db yu s i n gr a t i o - m e t h o dw er e m o v e dt h ea b n o r m a l e f f e c t so ft h et o p o g r a p h y ，a n dg o tt h er e s i s t i v i t yv a l u ew i t h o u tt o p o g r a p h i ce f f e c t s c a l c u l a t i n gs e v e r am o d e l sa n ds h o wt h er e s u l t st op r o v et h a tt h em e t h o do ff e m s i m u l a t i o na n dt h et e r r a i nc o r r e c t i o ni sa c c u r a t ea n df e a s i b l e k e yw o r d s ：f o r w a r d ；f i n i t ee l e m e n tm e t h o d ；n u m e r i c a ls i m u l a t i o n ；d i r e c tc u r r e n t r e s i s t i v i t ys o u n d i n g ；t e r r a i nc o r r e c t i o n i i 目录 摘要i a b s t r a c t i i 第l 章绪论1 1 1 概述1 1 2 国内外研究现状及问题的提出2 1 3 本论文的主要工作4 第2 章基本原理5 2 1 点源二维电场的特征5 2 2 边值问题8 2 3 变分问题1 3 2 4 有限单元法1 4 2 5 视电阻率的计算1 9 第3 章正演模型的计算2 3 3 1 水平层状大地的计算2 3 3 2 几个层状模型的对比计算及拟合情况2 5 3 3 复杂结构地电模型的计算2 8 第4 章地形改正3 4 4 1 模拟纯地形视电阻率计算3 4 4 2 模型算例3 7 第5 章结论与建议4 1 致谢4 2 参考文献4 3 附表4 5 个人简历、申请学位期间的研究成果及发表的学术论文4 8 桂林理工大学硕士学位论文 第1 章绪论 1 1 概述 电阻率法是以地下不同岩石和矿石的导电性差异为物性基础，通过观测与研究 人工建立的地中稳定电流场的分伟规律，以找矿石和解决各种地质问题的一组电法 勘探的分支方法。随着观测方式的不断改进以及解释水平的不断提高，其应用领域 在不断地拓宽。实践表明，电阻率法在普查金属、非金属矿产和研究地质构造方面， 以及在水文、工程、环境地质勘察等方面，都已取得了良好的地质效果，发挥其重 要作用。 一般地，为取得良好的地质效果，在电阻率法勘探中，常常根据不同地质任务 和不同地电条件，采用不同的装置类型，而装置类型是指一定的电极排列形式，由 于电极移动方式的不同，在电阻率法中又有电阻率剖面法和电阻率测深法之分。电 剖面法是采用固定电极距的电极排列，沿测线逐点地供电和测量，获得视电阻率剖 面曲线，以了解地下一定勘探深度内沿测线水平方向上岩石的电性变化，即用来探 查地下一定深度范围内的横向电性变化，以此解决多种地质问题；电阻率测深法的 全称为“垂向电阻率测深法”，顾名思义，它是用逐步改变供电电极距大小的办法 来控制勘探深度，由浅入深，以了解一个测点地下介质电阻率的垂向变化，可以综 合每条测线的各测点测量结果，通过定性和定量解释，即可获得地下岩石沿水平和 垂直方向变化的综合资料。 在电法勘探中电阻率测深是一种重要的方法，常用来解决不同电阻率地层的分 布问题，可用于探查油气、煤f f l 的地质构造，以及水文地质和工程地质的调查中， 其装置简单、操作简便、曲线易于解释等特点在我国得到广泛应用，在矿产、水资 源和工程勘查方面都发挥了巨大作用。电测深的装置有：通过改变a m 距离的二极 测深装置和三极装置；改变a m 和n b 间距离的联合三极测深装置；改变a b 间距离 的对称四极测深装置；以及改变偶极间隔常数的偶极测深装置。在实际的工作中， 一般采用对称四极测深装置，在施工条件限制时，也可采用三级测深装置，其他的 装置则比较少用。对称四极测深装置见图卜1 ： e 朋厂 h b 芒土 图1 1 对称四极测深排列示意图 桂林理工大学硕士学位论文 a b 为供电电极，m n 为测量电极，它们都对称于装置中心o 点，地面的测点和装 置的o 点重合。选择供电极距时，要求最小的极距应能反映第一层电阻率，最大的 极距能满足勘探深度的要求，并保证测深曲线尾支的完整，可解释最后一个电性层。 11 当a b 增大时，m n 也需要适当加大，通常要求m n 满足条件：三彳曰m n a b 【l 】 3 0j 测深剖面布置要尽可能的垂直地层走向，这样就能得到地下电性沿剖面和深度变化 情况。 在地球物理的正反演理论中，讵演模拟是问题的基础所在，地球物理反演是最 终的目的，但反演大多以j 下演模拟为基础。不同的正演模拟的手段，有着不同的速 度与精度，而且正演计算决定反演计算的速度、精度及其适应性和稳定性，因此， 所以选择一个兼顾反演稳定性和精确性要求的正演计算方法也是顺利实现反演的 前提。 1 2 电阻率法正演模拟的发展与研究现状 已知地电模型和场源分布，求解场的分靠规律，称为电法勘探的正演问题。正 演是电法反演问题的计算和观测资料解释的重要基础。一般求解电阻率法正演问题 有解析法、数值模拟法和物理模型法三种途径。 ( 1 ) 解析法，对地下形体比较规则的地质体，如球体和板状体、水平层等，可以通 过解析法推导出场值的解析表达式，该方法得到的结果很精确，但公式的推导过程 非常复杂，应用范围不宽，对于一些比较复杂的场源分布或不规则地质体，解析解 难以得到，不能满足实际生产的需要，一般用来验证其它正演方法的准确性。 ( 2 ) 物理模型法，主要使用各种物理模拟设备如水槽、土槽、电阻网络、导电纸和 薄水层等方法，它们均以物理现象的相似性为基础。在模拟无限大的大地时，模型 尺寸要比野外实物小得多，以同样的比例尺进行模拟，所用电阻率也不一定与野外 实际情况一样。只要保持地电模型按比例缩小，各地电体的电阻率比值保持不变， 便可用模拟试验得到与野外相似的电阻率曲线，这便是相似性原理。测量各种模型 的场值分布与大小，便可得到正演模拟结果。缺点和解析法一样，由于难以拟合比 较复杂的模型，因此应用范围也不是很广。 ( 3 ) 数值模拟法，根据地球物理中的偏微分方程和边界条件，用数值方法求解场值 的近似解，虽然是一种近似的方法，但由于它适用于复杂物性分布和复杂边界形状 的地球物理计算，所以适用范围很广。 电阻率二维正演数值模拟的方法中，最常用的大致分为三种：有限差分法、边界 单元法和有限单元法。 桂林理工大学硕士学位论文 有限差分法是一种经典的数值模拟计算方法，其基本原理是用差商代替微商， 把定解问题转化为代数方程组的求解。边界单元法的前身是边界积分方程法，随着 有限单元法的兴起，其单元划分和插值函数的概念引入到了边界积分方程法中，发 展成为边界单元法。有限单元法是将要分析的连续场分割为很多较小的区域，它们 的集合代表原来的场，然后建立每个单元上待求场量的近似式，再结合起来进而求 得连续场的解。 三种数值模拟方法在电法j 下演计算方面都有一定的优势和不足。有限差分法的 优点是方法简便易算，缺点是当物性参数复杂分布或场域的几何特征不规则时，适 应性比较差。而边界单元法的优势是j 下演速度快，内存需求少，主要用于地形改j 下 和地下少量地质体的正演模拟。有限单元法与前述方法相比，由于在处理复杂的几 何形状时，其灵活性和适应性比其它方法好，并且其计算步骤具备规范的流程，因 此在电磁法二维正演模拟中应用较多，特别是针对复杂地电断面的电磁场模拟其效 果更为突出。 1 9 7 1 年，j h c o g g o n 2 】首先将有限单元法引入到电法勘探领域正演模拟中，他从 电磁场总能量最小原理出发，实现了二维地电断面有限单元法正演计算，不过由于 有限元网格缺乏通用性，计算精度和速度未能达到实用水平；r o d i t 3 1 ( 1 9 7 6 ) 发展了 有限单元法的剖分方法，采用矩形网格剖分，以解决二维大地电磁测深j 下演问题， 使有限元向前发展了一大步；r i j o 4 1 ( 1 9 7 7 ) 引入了一种通用性网格剖分方法，使有限 元j 下演的精度和速度得到大幅度提高，成为计算二维地电条件下电法电磁法j 下演模 拟的有效工具，使有限单元法正式进入实用阶段。7 0 年代末，有限元法在甚低频法 正j ! 寅( e k a i k k o n e n t 5 1 ，1 9 7 9 ) 、时问域电磁法的正演模拟( j o n e t k u oe ta 1 【6 j ，1 9 8 0 ) d ? 得到 应用；d e ya n dm o r r i s o n 7 ( 1 9 7 9 ) 将混合边界条件引入到有限单元法的j 下演模拟中使 其得到进一步发展；1 9 8 1 年，p r i d m o r e 等【8 】将有限单元法用于解决三维地电断面的 电法电磁法正演问题。我国对有限单元法的引入( 朱伯芳p 】，1 9 7 9 ) 始于7 0 年代术。 8 0 年代初，我国数学家李大潜发表专著有限元素法在电法测井中的应用【l ，有 限单元法正式应用到电法领域，此后，地球物理学家周熙襄 i q t l 2 ( 1 9 8 0 ，1 9 8 6 ) ，罗延钟 【1 3 】( 1 9 8 7 ) 和徐世浙【1 4 】【1 5 】【1 6 】【1 7 】【1 8 】【1 9 】【2 0 ( 1 9 8 2 ，1 9 8 4 ，1 9 8 5 ，1 9 8 6 ，1 9 8 8 ，1 9 9 4 ) 等对有限单元 法进行深入研究，发表了一系列有限单元法在电法勘探中应用的论文或专著，研究 领域涉及直流电法、电磁法领域，网格剖分也由原来的简单剖分发展到三角单元和 三角一矩形综合剖分，周熙襄等在二维电阻率有限元法正演中采用混合边界条件等 优化措施，提高了计算精度和速度，对有限单元法的发展起到了推动作用；二维地 电构造中点源场的迭代有限元法( 杨进【2 1 1 ，1 9 9 3 ) 在点源二维电法j 下演的有限单元法 桂林理工大学硕士学位论文 ( 周熙襄【2 2 】，1 9 8 3 ) 一文的基础上从区域划分、网格剖分方面作了改进，将模拟区划 分为内区和外区，内区沿用矩形一三角形网格，外区则采用梯形一三角形网格；迭 代有限元法数值模拟若干问题的讨论( 杨进【2 3 1 ，1 9 9 4 ) 对上文的边界条件、网格剖分 以及模拟精度等问题作了分析：电导率分块线性变化二维地电断面电阻率测深有限 元数值模拟( 阮百尧【2 4 】，1 9 9 8 ) 一文则在总结前人工作经验的基础上又有创新，在单 元内，电导率、电位均采用双线性插值，计算结果精度较高。在有限单元i f 演模拟 方法的计算中，电源点是一个奇点，在计算中以有界函数来逼近，非常勉强地假设 在点源附近单元中v 是空间坐标的线性函数，所以点源附近的计算必然会存在误 差。且在有限元计算中，用傅立叶变换将三维问题转化为波数域的二维问题时，有 限的波数个数也会给有限单元法的计算带来误差。为此，解决这一问题将是本文的 一个重点。 1 3 本论文的主要工作 本文针对已有电阻率正演算法所存在的某些不足，结合当日i 一些较新的理论研 究成果，在前人研究的基础之上，作了一点探索性研究工作。本文工作主要包括以 下几个方面： l 、查阅国内外有关本课题的研究现状、分析了存在的问题与不足；构架了本 文工作的研究思路： 2 、结合已有的j 下演算法的理论与实践，进行点源二维电场的相关公式推导计 算，用有限单元法正演模拟计算； 3 、在5 个波数的有限元点源二维直流电测深总电位法模拟计算中，对电位进 行修正，计算几个相关模型，验证该方法的有效性。 4 、在有限单元法二维正演基础上，对带地形二维正演结果进行地形校正。 5 、对几个水平和起伏地形下的复杂结构模型进行计算，验证程序的正确可行 性。 4 第2 章基本原理 本章在点源二维电场的模拟的基础上，应用有限单元法进行正演模拟计算，在 j 下演过程中采用三角单元剖分网格，有利于对复杂地形的模拟，从而往后实测数据 在反演时不需要进行地形改j 下；电导率( 和电位在相邻节点间线性变化，使得电导 率和电位与实际地质构造更接近、与实际地电场变化规律更易吻合；采用双网格系 统，用大网格反映地下电性变化和测深点的位置，小网格进行实际的有限元数值计 算。 2 1 点源二维电场的特征 由电学方面的知识我们可以知道，表征物体导电性的物理量是它们的电阻率。 物体的电阻率愈小，它的导电性就愈好。根据实验结果，当电流通过柱状导体( 图 2 1 所示) 时，有如下关系式： 夕蚺孚 ( 2 1 ) 图2 - 1电流通过柱状导体体示意图 埋藏在地下的岩石和矿石，其电阻率是不能用( 2 1 ) 的关系加以测定的，但 是各种岩石和矿石之间的电阻率的差异可以在人工建立的电场中反映出来。为此， 我们将两个( 或两组) 金属电极a 和b 打到地下，使其上端分别与直流电源的正 极和负极相连。于是电流由其中一个电极输入地下，通过地下物质的传导，再由另 外一个电极流出地面，构成闭合回路，从而建立稳定的电流场。电极a 和b 称为 供电电极，一般情况下都可以认为它们是在供电时都是点电源。 2 1 1 单个点电源的电场i l j 如果电极距a b 非常大，那么在研究a 电源附近的电场分布时，b 电源产生的 影响可以忽略不计。这时可以认为b 电源位于“无穷远”处，因而由a 电源建立 的电场就是单个点电源的电场。在上述地面平坦、地下物质导电性均匀的条件下， 从a 电极进入地下的电流将由于空气是绝缘物质而全部在地下半空间内均匀地向 各方流动。 桂林理工大学硕士学位论文 我们可以用电流线表示电流分布。电流线上每一点的切线方向都与电流方向符 合；电流线愈密集，表示电流密度( 垂直于电流线的单位面积上通过的电流) 愈大， 电场愈强。在这个点电源的电场中，电流线应该是一系列由点电源出发，向各个方 向呈辐射状分布的直线，如图2 2 所示。 l i m i - 钇刚乒。l ( a ) 一一一 电流线 ( b ) 一一一等电位线 图2 - 2 电源点附近的电位和电流线示意图 另一方面，与电流线处处正交的面是电位处处相等的面，称为等电位面或等位 面，在这里应该是一系列以a 为中心的半球面( 在剖面图和平面图中只能分别表示 为一系列的半圆形曲线和同心圆，即上图中的实线) 。假定由a 点进入地下的电流 强度为i ，地下半空问的物质电阻率为p ，那么，对于地面任意一点m 而言，从理 论上可以写成该点的电位以为： u 0 ：j 丝 ( 2 1 1 ) 。 2 7 a m 式中彳肘为m 点与a 点间的距离。 根据( 2 1 1 ) 式，很容易求出地面上任意两点m 和n 之间的电位差为： u 知= u 二叫= 鲁( 去一六) 晓m ， 式中u ：是n 点电位，彳是n 点与点电源a 之间的距离。 6 桂林理工大学硕士学位论文 2 1 2 两个异性点电源的电场 在供电电极a 和b 各自建立的电场相互有影响而不能忽略的情况下，地下的电 场就应当是a 和b 各自引起的电场的叠加，也就是两个异性点电源的电场其分布 图2 3 所示： 嵝缈 、 一 ? ( a ) 剖面图 ( b ) 平面图 图2 - 3 两个异性点电源的电场模拟示意图 如果电流由a 电极进入地下而通过b 电极返回地上，就把a 当作是正电源( + i ) ， b 当作是负电源( - - i ) 。因此，根据( 2 1 2 ) 式及场的叠加原理可以得出，在两个 异性点电源的电场中，地面上任意两点m 和n 之间的电位差为： 噶= + = 篆( 菰1 一菰1 一丽1 + 丽1 ( 2 ) 式中丽、一a n 、丽、一b n 分别是a 与m 和n 以及b 与m 和n 之间的距离。令： 肚t t 竺丌 ( 2 - 14 ) ：= = = = = 一：= = = = ：= = = = = + ：= = = = = a ma nb m b m 并将嘴简化为【，可得：p ：k 华 ( 2 1 5 ) k 是仅与a 、b 、m 和n 四点的相对位置有关的量，称为装置系数或排列系数。 在单个点电源的情况下( 2 1 5 ) 式仍然适用，只是： k ：塾：2 万a m a n ( 2 1 6 ) 11 m n a ma n 桂林理工大学硕士学位论文 由( 2 1 5 ) 式可以知道，只要采用两个( 或两组) 供电电极a 和b ，使地下建立 人工的稳定电流场，已知电流强度为i ，并用两个测量电极在平坦的地面上任意两 点m 和n 打入地下以观测电位差u ，就能在地下物质电性均匀的条件下求出该物 质的电阻率p 。 2 2 边值问题 2 5 】 在地面彳点，置一个电流强度为i 的点电源，用表示电流密度矢量法。在空 间作任意闭合面厂，q 是厂所围的区域。根据通量定律，当彳点在闭合面厂上时( 如 图2 4 a ) ，流过闭合面厂的电流总量为，；当彳在闭合面厂以外( 如图2 4 b ) ，流 过闭合面厂的电流总量为零，写成数学式子： 妒= 仨筹 2 上面是面积分，为简单起见，此处的面积分用单积分号表示。 k 2 泌 1 4 田 y 锹 - _ _ ， 弋 图2 - 4 模型边界示意图 ( b ) 根据数学分析中的奥高公式，可将矢量的面积分转化成矢量的散度的体积分： = 百v dq = ：三：? ql 一 以j 翻夕表示以a 为中心的每函数，根据j 函数的积分性质，有： 呵万ca ： 王乞： q l2 将前面两式的对比，得 v j = 2i 万( a ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) 桂林理工大学硕士学位论文 电流密度，与电位梯度有如下关系： j = 一仃vu ( 2 2 5 ) 其中6 是介质的电导率，将( 2 2 5 ) 式带入( 2 2 4 ) 式，得电位应满足的微分方程： v ( 仃vu ) = 一2 ，艿( 彳) ( 2 2 6 ) 在直角坐标系中，将上式展开，得 昙( 盯孚) + 晏( 仃罢) + 拿( 盯孚) ：一2 1 6 ( x a ) 万( y 爿) 万( 乃) ( 2 2 7 ) o xo x卯o zo z 其中x a ，y a ，z a 是彳点的坐标。( 2 2 7 ) 式是三维微分方程，对于二维构造，选取如 图2 5 所示坐标系，使z 轴平行构造走向。由于6 在z 方向不变，昙p 罢) ：仃鲁， 但是否在x 、y 方向变化，6 不能提到微分符号旦旦之外，( 2 2 7 ) 式变为： 苏砂 拿( 盯罢) + 昙( 盯罢) + 拿( 仃罢) ：一2 1 6 ( x a ) 万( j ，) 万( 乙) ( 2 2 8 ) o x呶 0 、)o v o zo z 这是二维构造中点电源电场的电位所应满足的微分方程，也是三维的点电源电场微 分方程。 x 图2 - 5边界示意图 现用一个足够大的区域代替半无限空间，在地下作一个半径足够大的圆柱面 r 。，它与地面f 。组成一发封闭的区域q ，r 。+ r 。是q 的外边界，外边界条件是： 1 ) 、在地面r 。上，电流沿地表流动，所以电流密度的法向分量为零，也就是电流 的法向导数为零 o ul ，：0 o n 1 ， 9 ( 2 2 9 ) 桂林理工大学硕士学位论文 2 ) 、在无穷边界r 上，假定区域q 内部的电性不均匀性对凡上电位不发生影响， 其电位是点电源电位，即 i c “k2 7 ( 2 2 1 0 ) 其中c 是比例系数r 是电源a 至边界1 1 。上的距离，( 2 2 1 0 ) 式中的假定是近似的， 一般来讲，愈大，近似性愈好，所以在用有限元解无限区域的问题时，需要将区 域取的足够大，( 2 2 1 0 ) 式是第一类边界条件，其中含有未知常数c ，下面将会看 到，将它变成第三类边界条件，可以消去常数c 。 在区域内部，存在着电性不同的介质的分界面，这些边界称为内边界。内边界 条件是自然边界条件，在泛函极值过程中，将会自动得到满足，在点源二维电场中， 可以证明内边界条件也是自然边界条件，后面的讨论将不再涉及内边界条件。 现在将二维构造中点电源电场中的电位的边值问题( 简称为点源二维问题) 归 纳如下： 拿( 仃孚) + 晏( 盯豢) + 盯石0 2 u ：一2 1 8 ( _ ) 万( 儿) 艿( 列q o xo x c c 、 o z 孰_ 0 l ( 2 2 1 1 ) “卜争叽 这是三维问题，对于这种类型的三维问题，可以用傅立叶变化将它变成二维问 题，用“伍只矽表示三维空间中的电位，对u 在z 方向进行傅立叶变换： u ( x ，y ，k ) = l u ( x ，y ，z ) c o s k z d z ( 2 2 1 2 ) 5 简记作泸厅训，因为在点源二维问题中，电位“阮只矽对于x , y 平面，即“肜矽= “伍彭- z ) ，所以( 2 2 1 2 ) 式中用余弦进行变换，积分区问从0 到。 1 ) 微分方程的变换 根据傅立叶变换的性质，有： f 旦仃塑】一旦仃塑f ；仃丝】- 旦莎型 la x a x a x8 x l 卸卸。巩如 f i 害i 一2 u 1 0 桂林理工大学硕士学位论文 根据。函数的积分性质，有 了。) 万( y a ) 6 ( z 。：导6 ( _ ) 艿( y a ) ：了1 8 ( x ) c o s k z d z 万( 彳) i 。) 万(6 ( z 。 = 寺6 ( ) 艿( y a ) = i 万( 彳) 式中，a 表示a 点在x y 平面上的坐标，对( 2 2 1 1 ) 式的微分方程的两端进行傅立 叶变换，得： 去( 仃警) + 导卜詈卜后2 训= 埘c 彳， c 2 2 t 3 ， 这个方程只含有x 、y 坐标，所以它是二维微分方程，但是方程中含有瓜它起到参 数作用。用二维哈密顿算子v = 晏巳+ 昙p ，将( 2 2 1 3 ) 式简记为 v ( 仃v u ) - k 2 t r u = - ：8 ( a ) 前面三维空间中所用到的v 是三维的，现在已经将三维问题转化为二维问题，和下 面所用的v 都是二维的。当然这只是单点电源的情况，双点电源的话只需在( 2 2 1 4 ) 右边源的地方加上b 点电源的舾) 函数，即为： v ( 口v u ) - k 2 c r u = 一，( 万( 彳) + 万( b ) ) q( 2 2 1 5 ) 2 ) 边界条件的变换 ( 1 ) 要l = o 的变换 法向n 处在x y 平面内，而且与z 无关，所以 榭崩r j _ 。 2 舶， ( 2 ) “i r 号一号的变换 u = f 昙一丢】2 ，了每一了季虿，c o s k z d z = c k o ( 饥) - 丛。( ) ( 2 2 1 7 ) 桂林理工大学硕士学位论文 其中r 是通过坐标原点与走向垂直的截面内，f ，边界上的点至坐标原点的距离，k o 是第二类零阶修正贝塞尔函数，k o ( x ) 的微商是 旦d x k 。( x ) = 一k ( 工) 其中k 。是第二类一阶修止贝塞尔函数。k o ( x ) ，k i ( x ) 的图形见示意图2 - 6 0 123 45 图2 - 6 第二类零阶、一阶贝修正塞尔函数曲线 对于( 2 2 1 7 ) 在法向1 1 方向进行微商，得 _ o u ：c _ o k o ( k g ) i o r a o k o 。( k r b ) o 刍r b 一二涨i ( 矾) c 。s ( 么，玎) + c k k 。( 饥) c 。s ( ，胛) 锄 西。 锄 钆 锄 ”“ 1。 ( 2 2 1 8 ) 其中c o s 化，z ) 是矢径，与外法向靠的夹角的余弦，从( 2 2 1 7 ) 式解出c ，代入上式， 消去c 后，得 降o n 七丝巡k o 跨铲k ou 卜 2 加， l( 帆) 一( )j i r 上式是第三类边界条件。 综上所述，经过傅立叶变换后，三维边值问题( 2 2 1 1 ) 变换成t y l j - 维边值问题 1 2 5 4 3 2 1 0 桂林理工大学硕士学位论文 v ( o v v ) - k 2 仃u = - i a r ( a ) + i a ( b ) q 型：o r o n 5 _ o u + kk , ( k g ) c o s ( g , n ) - k l ( k r e ) c o s ( g , n ) u ：0 r j l 。 一- r t 锄 k o ( 帆) 一k o ( 饥) 。 ( 2 2 2 0 ) 前面，( 2 2 i i ) 式子中的区域q 是三维的，边界r ，f 。是三维区域的边界，此处， ( 2 2 2 0 ) 式的区域是二维的，边界1 1 ，f 。是二维区域的边界。 2 3 变分问题 2 5 1 用有限单兀法求群边值i 司越( 2 2 2 0 ) ，首先妥将它变成燹分i 司题。f 面，我们 给出于这个边值问题相应的变分问题。 构造一个泛函： 州) 。【詈( 跗) 2 + 圭旷彬小) u + i o - ( 印 q ( 2 3 1 ) 其中积分区域q 由地面边界r ，及无穷远边界r 。围成。其变分 o f f ( u ) = v u v 6 u + k 2 0 u d u 1 6 ( a ) d u + 1 6 ( b ) d u d n = i v ( o v u d u ) d q + 一v ( 田u ) + 尼2 盯u 1 6 ( a ) + 1 6 ( b ) 6 u d n 根据微分方程( 2 2 2 0 ) 可以知道，上式右侧最后一项的被积函数必为零，所以 刖= q v ( o r6 u ) d f ：= ，哩盯珈r ( 2 3 2 ) q r 。+ r 。 式中的n 代表q 的外法向，将( 2 2 2 0 ) 式中的外部边界条件代入( 2 3 2 ) 式，得 万，( u ) 2r 仃o - - 筹d u d r = 一矾! 量鱼壁2 蔓里薹善弓言 三鬟杀铲u 6 u 万r ：一万i ! 陬塑蚴竺虹尘墅业竺幽u z d fl i 2 r j 【0 ( 饥) 一k o ( 饥) l 或 仃l ，c u ，+ 丢r j 仃足鱼l 垡兰上兰善老吉弓 j 詈轰笔产u2 d r 1 桂林理工大学硕士学位论文 = 万 詈c v ，2 + 丢k 2 0 u 2 - 1 6 ( 彳) + ，万( b ) d q + 上f 盯后墅匦也业也立坐堕型亟盟uz d f ：0 2l k o ( 饥) 一k 。( k r b ) 所以二维边值问题( 2 2 2 0 ) 与下列变分问题等价 f ( u ) = 詈( v u ) 2 + 圭k 2 0 u 2 - 1 6 ( 彳) ? + ，万( b ) u p q + ! f c r kk , ( 帆) c o s ( r a , n ) - k , ( k r b ) c o s ( t o , n ) u 2 d f ， 2 一 k o ( 饥) 一k o ( 蛾) 6 f ( u ) = 0 2 4 有限单元法 1 】 ( 2 3 3 ) 将积分区域剖分成许多三角单元，总节点数为n ( 图2 7 ) ，则( 2 3 3 ) 式中对 区域q 的积分可分为对各三角单元的积分。 f ( 【，) = 丢萋 盯 ( v u ) 2 + k 2 u 2 一2 ，万( 彳) u + 2 1 6 ( b ) u ) d + 三f 盯尼益墨訾要掣缨， q 4 j 2t x o ( 峨) 一k o ( 饥) ” 2 4 1 线性插值 图2 7 区域q 有限元剖分示意图 图2 8 ( a ) 是母单元，( b ) 是子单元，两个单元的图示如下所示： 1 4 桂林理工大学硕士学位论文 o ( a ) 三角形母单元( b ) 三角形子单元 图2 - 8 三角单元设计图 设电位u 和电导率盯是单元内的线性函数： u = n , u l + n 2 2 + n , u 3 0 - = 1 0 - t + n 2 0 2 + n 3 0 - 3i ( x 。，y 。) x e 9n 1 、n 2 、m 称为形函数： l = 去( 叩+ 包y + q ) 2 = 去( 彬+ 6 2 y + c 2 ) 3 = 去( 叩+ 岛y + c 3 ) ( 2 4 2 ) 口l = y j 一虼6 i = 一x jc i = x y m x y j a 22y m y ib 2 = x i x m c 2 = x m y i x i y ” a 32y t y j 6 3 = x j x ic 3 = x i y j x j y i = 三( q 包一以：6 1 ) 2 4 2 单元分析 ( 2 4 3 ) 将( 2 3 1 ) 中的区域积分，分解为各单元e 的积分，积分的第一部分： 蔓巾u ) 2 d r 2 = 1 u r ek , 。玑泣4 舢 其中的- l 。是，1 1 0 - ( 、7 “) 2c ，( 2 ：= ；f ：i c r ( 罢) 2 - ( 考) 2 = 扛( 乃) 铲1 k 。， 将2 4 2 式代入得： 1 5 桂林理工大学硕士学位论文 k 。= m q + 弘巳+ 虬 8 n a x a n ， i g x a n d x i 盟盟盟| + l 玉 缸缸l 化简该积分后就得到： k 。= m q + 以巳+ 出咖 【q + 盯，+ ) = - - - - - - - - - - - - - 二- - - - - - - - - - - 一 1 2 q 2 q i 2 a m 2 a n