（地球探测与信息技术专业论文）基于ansys的激发极化法中地形影响问题的数值模拟.pdf

中南大学硕士学位论文摘要 摘要 目前，电法勘探在我国国民经济和社会发展各个领域应用非常普 遍，遇到的问题也多，特别是地形影响问题。我们知道，虽然地形不 平对激电异常的影响比对导电性异常的影响小，但为了正确推断解释 激电异常，仍必须考虑地形影响问题。目前地形改正的方法很多，本 文应用a n s y s 有限元进行大量的电阻率法和激电法数值模拟，然后用 比值法进行地形改正。众所周知，有限元法适合于模拟形状不规则的 模型，计算精度高，是地球物理数值模拟中较常用的方法。 本文首先回顾了电法勘探地形影响问题的研究现状，总结了现有 激发极化法的数值模拟方法，简单介绍了a n s y s 有限元软件的主要功 能，在进行激发极化法数值模拟过程中的主要步骤及其注意的关键问 题。然后介绍了点源场电位边值、变分问题；通过使用a n s y s 本身提 供的a p d l 语言编制相应的计算程序，以山脊、山谷地形为例，进行 了大量的电阻率法和激发极化法地电模型计算，并对模拟计算出来的 结果加以验证和分析，总结了起伏地形情况下的异常分布规律，得到 了一些比较满意的结果。 最后总结全文，表明了a n s y s 有限元软件在地球物理领域里有广 泛的发展前景。 关键字：激发极化法，有限元，地形影响，数值模拟 中南大学硕士学位论文 a b s t r a c t a b s t r a c t c u r r e n t l y ，e l e c t r i c a le x p l o r a t i o ni sv e r yw i d e l ya p p l i e di nv a r i o u s f i e l d so fn a t i o n a le c o n o m ya n ds o c i a ld e v e l o p m e n t h o w e v e r , t h e r ea r e m a n yp r o b l e m se n c o u n t e r e d e s p e c i a l l yt h et e r r a i ne f f e c t 。w ek n o wt h a t t h et e r r a i ne f f e c to ni pa n o m a l yi sm u c hl e s st h a nt h a t0 1 1c o n d u c t i v i t y a n o m a l y h o w e v e r , t h et e r r a i ne f f e c tm u s tb ec o n s i d e r e di no r d e rt o e s t i m a t ea n di n t e r p r e t1 1 a n o m a l yc o r r e c t l y p r e s e n t l y ，t h e r ea r em a n y m e t h o d so ft e r r a i nc o r r e c t i o n am a s so fn u m e r i c a ls i m u l a t i o n sf o r r e s i s t i v i t ya n di pm e t h o d sa r em a d eb yu s i n ga n s y si l lt 1 1 i sa r t i c l e w h i c hi sas o f t w a r ef o rf i n i t ee l e m e n ta n a l y s i s t h e ni ti si n t r o d u c e do f t h et e r r a i nc o r r e c t i o nb yr a t i om e t h o d a sw ek n o w t h ef i n i t ee l e m e n t m e t h o ds t l i t st os i m u l a t em o d e l so fi r r e g u l a rs h a p ea c c u r a t e l y i ti st h e c o m m o n l yu s e dm e t h o di ng e o p h y s i c a lh u m e r i c a ls i m u l a t i o n i nt h i sp a p e r , w ef i r s t l yr e v i e wt h es t u d ys t a t u so ft e r r a i ne f f e c ti n e l e c t r i c a le x p l o r a t i o na n ds u mu pt h e e x i s t i n gn u m e r i c a ls i m u l a t i o n m e t h o d so f 碑a n dt h e nw ei n t r o d u c et h em a i nf u n c t i o n so ft h ef i n i t e e l e m e n ta n a l y s i ss o r w a r e a n s y s 。t h ek e ys t e p sa n di s s u e si ni p n u m e r i c a ls i m u l a t i o n a f t e r w a r d ，w em a k eas i m p l ei n t r o d u c t i o na b o u t t h eb o u n d a r yv a l u ea n dv a r i a t i o np r o b l e m so fp o i n t s o u r c ef i e l d m a n y g e o e l e c t r i cm o d e l so fr e s i s t i v i t ya n di pm e t h o d sa r ec a l c u l a t e db yu s i n g t h e c o r r e s p o n d i n gp r o g r a m m e s i na p d ll a n g u a g eo f r e r e db yt h e s o f t w a r ea n s y si t s e l f , t a k i n gt h et e r r a i n so fr i d g ea n dv a l l e ya st h e e x a m p l e s 耶1 er u l eo fa n o m a l yd i s t r i b u t i o no nt h ec o n d i t i o no ff l u c t u a t e d t e r r a i ni ss u m m a r i z e db yv a l i d a t i n ga n da n a l y z i n gt h ec a l c u l a t e dr e s u l t s a 1 1a b o v e 。w eo b t a i ns o m es a t i s f a c t o r yc o n c l u s i o n s a tt h ee n d ，w eh a v ep r e s e n t e dt h es u g g e s t i o nt h a ta n s y sc a nb e w i d e l yu s e di ng e o p h y s i c a lf o r w a r dm o d e l i n gw i t hv a l u a b l ea p p l i c a t i o n s k e yw o r d s ：i pm e t h o d ；i n f i n i t ee l e m e n t ；t e r r a i ne f f e c t ；n u m e r i c a l s i m u l a t i o n i i 原创性声明 本人声明，所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。尽我所知，除了论文中特别加以标注和致谢的地方外，论文中不包含 其他人已经发表或撰写过的研究成果，也不包含为获得中南大学或其他单位的学 位或证书而使用过的材料。与我共同工作的同志对本研究所作的贡献均已在论文 中作了明确的说明。 作者签名 日期：年一月二日 关于学位论文使用授权说明 本人了解中南大学有关保留、使用学位论文的规定，即：学校有权保留学位 论文，允许学位论文被查阅和借阅；学校可以公布学位论文的全部或部分内容， 可以采用复印、缩印或其它手段保存学位论文；学校可根据国家或湖南省有关部 门规定送交学位论文。 作者签名起沥导师签名：嫩 中南大学硕士学位论文第一章绪论 1 1 引言 第一章绪论 在地球物理的理论计算中，存在着两类基本问题嘲：正问题和反问题。 给定场源的分布，求解场值的大小，这是正问题，或称正演问题。例如，在 磁法勘探中，给定了铁矿体的分布形态、磁化率和磁化强度，计算铁矿体产生的 磁异常；在电法勘探中，给定了各种岩石电阻率及岩石的分布形态，计算与各种 测量排列装置相应的视电阻率；等等。 长期以来，电法正演问题的计算中用场论中的经典方法解析求解，都只针对 于少数简单规则的地电条件。除此以外，对大量的实际地电条件尚未得出解析求 解的公式。 解决地球物理正问题的途径【7 l 有：( 1 ) 解析法，即用解析公式表达场值的大 小，这是各种地球物理教科书讲的最多的方法。但是，只对少数规则形体，例如 球体，板状体和水平层等，才能用解析法导出场值的解析表达式，因而，此方法 实用的范围非常有限。( 2 ) 模型实验法，这是电法勘探中使用较多的方法。例如， 各种水槽、土槽模型实验。但模型的制作比较麻烦，尤其是物性分布复杂的模型 很难制作。( 3 ) 数值模拟法，即根据地球物理中的偏微分方程和边界条件，用数 值方法求场值的近似解。这是一种近似方法，但它实用于复杂物性分布和复杂边 界形状的地球物理计算，所以适用范围非常广泛，数值模拟法需要进行大量的计 算工作。近三十年来，由于计算机技术的迅速发展，克服了计算工作量浩大的困 难，数值模拟法已成为地球物理正演的最主要的方法。 在大多数勘探运用中，由于都在较平地形上进行，故可以不考虑地形影响。然 而在地形起伏较大时，如铁路、公路隧道勘探中，它导致异常曲线畸变严重【13 。 有时淹没地质体的信息，有时歪曲有用信息，有时呈现假信息，使人们不能做出正 确的推断解释，降低了异常解释精度，甚至导致错误的解释。因而，研究地形影响 及其改正方法，是提高山区电法勘探地质效果的一个关键。本文通过对有地形影 响的异常体模型进行电阻率和极化率的数值模拟，并对所模拟的结果用比值法进 行地形改正，从中找到异常体的空间状态与曲线形态的内在规律性，对野外工作 及野外数据的定量解释具有一定的指导性。 1 2 电法勘探地形影响问题国内外研究现状 地形影响问题是2 0 世纪5 0 年代伴随着直流电法在矿产资源勘查等方面的大 规模应用开始受到重视。此后随着大地电磁法( 灯) 和可控源音频大地电磁法 中南大学硕士学位论文第一章绪论 ( c s a m t ) 为代表的无源和有源电磁法应用的扩展，人们也开始关注电磁法中的地 形影响问题。理论与实践表明，地形不但可以引起假异常，而且会掩盖地下由矿 体或目标物引起的真异常，如果不能正确的认识和消除，就会导致完全错误的解 释结果。 早期，日本、前苏联及我国，乃至西方国家的地球物理学者曾采用导电纸模 拟、电阻网络模拟、水槽模拟、土槽模拟等物理模拟方法研究地形影响问题，如 按一定比例制作导电纸、电阻网络、水槽、土槽等模型来模拟野外情况，得到地 形影响的定性认识，并在野外实际测量中加以识别。有些学者也应用模拟结果， 对野外观测数据作些简单的计算和改正，试验将地形影响去掉。这些做法基本属 定性解释。其后，也有学者采用一些解析计算方法研究地形影响，如保角变换法 和基于保角变换的坐标网法、角域地形叠加法、格林理论、模拟电荷法等。但由 于数学公式的复杂性和只能计算简单地形，使应用受到定限制。 近些年来，随着计算机技术的发展和物探解释方法的进步，国内外学者在水 平地形的电法、电磁法资料2 - d 、3 - d 数值模拟和解释技术基础上，研究开发出适 合于起伏地形的2 - d 、3 - d 数值模拟和解释技术，从而使地形问题的研究更加方便 可行，也使复杂的2 一d 、3 - d 地形影响问题的模拟、识别和校正成为可能。起初是 用比值法进行校正，后来发展到带地形进行2 - i ) 、3 书反演直接消除地形影响，取 得了一定效果，定程度上提高了电法勘探数据的解释水平。 国内外电阻率、激发极化法、充电法等电法地形影响问题的文献比较多，尤 在电阻率2 - d 和3 - 1 ) 地形模拟与解释技术方面取得较大进展。 物理模拟方法：日本t a k e s h i ( 1 9 5 2 ，1 9 5 3 ) 借用物理模拟方法研究电阻率曲线 的地形影响。前苏联v a e s h e v ( 1 9 5 9 ，1 9 6 0 ，1 9 6 5 ) 借助物理模拟研究电阻率测量的v 型地形效应。在2 0 世纪6 0 7 0 年代，我国徐世浙、杨文采、易永森、葛为中、何 继善等曾采用导电纸2 一d 模拟、电阻网络2 - 0 或3 - d 模拟、水槽2 - d 模拟、土槽3 - d 模拟等物理模拟方法研究地形影响问题。这一时期基本是定性解释阶段。 解析计算方法：e h a n t u r i s h v i l i ( 1 9 5 7 ) 从理论上得到棱形和阶梯形不规则地 形的均匀电场的畸变。s p i e g e l 等( 1 9 8 0 ) 使用s c h w a r z c h r i s t o f f e l 保角变换方法 模拟2 - d 地形下均匀大地中含3 一d 局部小空洞的地形变化。k o s t y a n e v ( 1 9 8 5 ) 和徐 世浙等( 1 9 8 8 ) 使用格林理论于电阻率地形模拟，徐世浙等( 1 9 8 8 ) 的研究与 k o s t y a n e v ( 1 9 8 5 ) 在结果方程上有所不同，k o s t y a n e v 直接使用面积分公式，而徐 世浙等使用立体角方法导出位方程。刘阳等( 1 9 9 6 ) 将模拟电荷法引人纯地形异常 的计算，计算几个典型的2 - d 角域纯地形异常及一个3 一d 坡顶的纯地形异常。 b i s w a s b h a t t a c h a r y a ( 1 9 9 8 ) 针对埋在均匀起伏地形里具有明显电阻率和极化 差异的2 - 1 ) 柱极化体提出由线性偏微分方程主导的2 加n e u m a n n d i r i c h l e t 型边 2 中南大学硕士学位论文第一章绪论 界值问题的一种精确解法，研究起伏地形下深度、隆起地形的宽高、电阻率比等 对电场振幅的影响，采用“比值法”进行地形校正。所谓比值法即是用实测视电 阻率除以“地形校正因子( 地形模型视电阻率地形模型电阻率) ”，从而得到改 正后或平地形下的视电阻率。 面积分方程法：g l op p l i g e r ( 1 9 8 4 ) 对充电法和磁电阻率法( 埘r ) 3 - d 地形影 响的模拟利用面积分方程方法实现。开发出面电荷源像的新方法处理空气一地球 界面电场边界条件，用有限个矩形面元模拟地形表面，面元各自倾向和电荷密度 不同。电荷值用迭代解技术得到。面电场根据面电荷分布、电流源和像来计算。 磁殇由涉及面倾向和电场的面积分计算。展示了充电法和删r 野外测量数据地形 校正实例。 边界元方法【：徐世浙( 1 9 8 4 ) 介绍起伏地形上2 5 - d 电场计算的边界元方 法。黄兰珍等( 1 9 8 6 ) 介绍点源场电阻率法2 一d 地形改正的边界元法，并计算了山 脊、山谷、角域地形上联合剖面法视电阻率异常曲线，展示用比值法对实测曲线 进行地改的应用效果。徐世浙和汪晓东( 1 9 8 7 ) 介绍用边界元法计算争d 山脊地形 下中梯法低阻球体视电阻率异常的方法和用比值法进行地改的算例。徐世浙等 ( 1 9 8 9 ) 对2 - d 电阻率地形改正边界元算法做了改进，减少了计算量，用比值法对 对称四极剖面曲线进行了地改。钱家栋和马钦忠( 1 9 9 2 ) 讨论用边界元法解决层状 介质中存在2 一d 横向不均匀情况下点源场视电阻率计算问题。徐世浙( 1 9 9 3 ) 介绍 边界元法模拟点源2 - 0 地形影响的结果。黄兰珍和田宪谟( 1 9 9 7 ) 介绍边界单元法 应用于工程地质勘查电阻率法的地形改正实例【4 8 】。刘继东( 1 9 9 8 ) 介绍求解异常 电位的2 - d 边界单元法对电测深资料进行比值法地形改正的应用效果。徐世浙等 ( 2 0 0 2 ) 用2 - d 边界元方法讨论2 _ d 地形对经向剖面的视电阻率影响。田宪漠等 ( 1 9 8 6 ) 介绍点源电阻率法3 一d 地形改正的边界元法，并计算了3 一d 山脊地形上联剖 的以n 值，对相应地形下埋藏着良导电性椭球体的p 。值用比值法进行了地形改 正。徐世浙等( 1 9 8 8 ) 提出模拟3 - 0 地形影响的新方法，利用格林函数将3 - d 边值问 题转换成积分方程问题，用边界元方法解积分方程：给出用比值法电阻率3 _ d 地形 校正的例子。 有限差分法：吴小平等( 2 0 0 1 ) 开展起伏地形条件下有限差分电阻率激发极 化3 - d 正、反演技术研究瞄】。 有限元方法t 2 0 , 4 6 ：c o g g o n ( 1 9 7 1 ) 针对直立导电极化体，用有限元法计算得到 平底、斜边的山谷和山脊2 _ d 地形下电阻率及极化率响应曲线。c o g g o n ( 1 9 7 3 ) 用有限元技术模拟比较极化率双偶极和三级装置的地形影响。r i j o ( 1 9 7 7 ) 使用有 限元方法模拟电和电磁数据。r i c h a r dcf o x 等( 1 9 8 0 ) 利用2 一d 有限元算法做了源 于地形的偶极一偶极视电阻率异常和地形( 山谷、山脊) 对极化率影响的系统研究。 3 中南大学硕士学位论文 第一章绪论 对视电阻率地形影响校正使用有限元法计算校正系数，用比值法进行校正。周熙 襄等( 1 9 8 6 ) 应用有限元法计算偶极装置下台阶、山脊、山谷2 一d 地形对视电阻率 异常曲线的影响及在地表下存在低阻体时山脊、山谷2 一d 地形的组合视电阻率异 常曲线。p a n a g i o t i sit s o u r l o s 等( 1 9 9 9 ) 用2 5 - d 有限元模拟方法考察几种常用 的电阻率装置的地形影响并采用比值法进行校正。杨华和李金铭( 1 9 9 9 ) 以山脊、 山谷地形为例，通过点源2 _ d 有限元数值模拟给出了起伏地形条件下地形对充电 法电场的影响，并用比值法对地形影响进行校正田j 。z o n g e 公司( 2 0 0 1 ) 推出有限 元2 - d 带地形商业反演软件( t s 2 d i p 、r4 1 0 ) ，可反演带地形电阻率r e 数据的电 阻率模型断面。ht r u m a nh o l e o m b e 和g e o r g erj i r a c e k ( 1 9 8 4 ) 给出3 一d 有限元 算法，对该算法进行多种测试并展示实际模型精确结果，采用比值法校正。 本文主要采用有限元法进行大量的激发极化发地形影响问题的数值模拟。 1 3 本文的主要工作 本文的工作主要由以下部分组成： 一、绪论：主要介绍了电法勘探地形影响问题的研究现状。 二、激发极化法及其数值模拟：首先简要介绍激发极化法的原理和进行 数值模拟的方法。然后对a n s y 8 软件及其平台进行简单的介绍，然后主要阐述总 结了在a n s y s 建模及其网格划分中注意的问题，应遵循的一些基本原则。 三、a n s y s 在有限元模拟中的实现：分别简要推导了点源三维电场、点源 二维电场问题的基本方程、边值问题及其变分问题，然后分析了如何应用a n s y s 软件实现激发极化法有限元模拟。 四、激发极化法有限元模拟与分析：对点电流源和线电流源激励的电场 进行了大量的有限元带地形的地电模型正演模拟。其中涉及到了有限元计算中的 源、边界条件、计算区域的选择、单元剖分等问题，应用a n s y s 编制a p d l 语言 程序，对激发极化法地形影响问题进行大量的数值模拟并加以分析。 五、总结全文：总结了应用a n s y s 实现激发极化法地形影响问题数值模 拟的一些结论和分析结果，说明a n s y s 在地球物理中必定会有广阔的应用前景。 4 中南大学硕士学位论文第2 章激发极化法及其数值模拟 ， 第2 章激发极化法及其数值模拟 2 1 激发极化法原理 激发极化法是5 0 年代和6 0 年代初，在我国开始实验研究和推广的。实践证明， 它是应用最广的和效果最好的一类电法勘探方法，早期是以直流( 时间域) 激电 法为主，通过长期的应用和研究取得了许多重要的成果。为了提高激电法的抗干 扰能力和减轻装备，7 0 年代初又开始研究推广了交流( 频率域) 激电法，这一时 期应用最多的主要是变频法。中南工业大学提出的双频道激电法是对变频法的发 展，为了异常的区分和电磁耦合问题，8 0 年代初又开始对频谱激电法进行研究， 取得了不少有价值的成果。 激发极化【2 】( i n d u c e dp o l a r i z a t i o n ，缩写为i p ) 是发生在地质介质中，因 外电流激发而引起介质内部出现电荷分离，产生一个附加的“过电位”( o v e r v o l t a g e ) 的一种物理化学现象。这个现象可用图2 1 的装置来观察。 图2 - 1 激发极化法现象示意图 f i 9 2 - 1t h es k e m h p l o t o f l p p h e n o m e n a 经由供电极a b 向地下供入直流电流i ，观察测量电极洲问的电位差及其随 时间的变化。如果地下介质是均匀非极化的，并且保持供电电流不变，那么在 心间测到的电位差形将如图2 2 中的水平虚线所示，也不会随时间变化。如果 地下存在可极化的介质，比如像图2 - 1 那样的金属矿体，则供入的电流流过矿体 时，矿体将产生电化学反应而被“极化”成原电池。此原电池也将产生电流，俗 称“二次电流”。二次电流在空间也要产生电场，在间形成电位差，称为二次 电位差。巧是随时间变化的，常记为( f ) 。供电时在删间所测到的电 位差a v ( t ) 是一次电位差a k 和二次电位差k ( ，) 之和，俗称总场电位差，它的 涵义是供电电场与极化电场二者在蜊问产生的总电位差。由于圪( n 随时间而 变，是时间的函数，故a v q ) 也是时间的函数，上述现象可用图2 - 2 来表示。 5 中南大学硕士学位论文 第2 章激发极化法及其数值模拟 图2 - 2 在时闻域观察激发极化法示意图 f i g2 - 2t h e s k e t c hp l o t 自f i pp h e n o m e n ai nt i m ed o m a i n 矿体( 或者别的可极化介质) 被极化是一个过程。过程随时间变化的速率 先快后慢，最后趋于饱和。饱和值的大小圪( s 是s a t u r a t e 的缩写) ，以及达 到饱和经历的时间t 与介质的性质有密切的关系，一般t 。为数秒到数分钟。 如果在图2 1 中把供电回路的开关k 断开，便会看到，虽然经由占供入地 下的电流已不复存在，但矿体的极化依旧存在，它将继续放电。因切断了激发电 流，随着放电的进行，矿体的极化也逐渐减小。相应地，断电后供电引起的电位 差k 已不再存在，但极化引起的电位差( ，) 并不立刻消失，仍然能观测到， 并且先快后慢地逐渐减小，经过数秒甚至到数分钟，才最终归于消失。如图2 2 中的曲线b 。上面观察到的现象有点像电容器的充放电，故习惯上将供电时心 间测得的总场电位差矿( r ) 的曲线称为充电曲线，而把断电后测到的二次电位差 巧( f ) 的曲线称为放电曲线。 上述电现象因外电流的激发而引起( 外电流接通以前决不会发生) ，外电流 切断后却并不立即消失，说明电流的激发导致了介质本身的极化，故称之为激发 极化。 2 2 激发极化法数值模拟 模型和模拟实验( 统称物理模拟) 方法，虽然比数学解析法更易于解决比较 复杂地电条件的问题，但仍有其缺点。比如，对比较复杂的地电条件构建模型也 十分复杂，此外，实验结果受电位差观测误差的影响，往往精度不高；最后，物 理模拟的效率较低，并且实验结果不便于用于自动化数据处理。现代数值模拟方 法可以避免这些缺点。 电场的数值模拟，就其原理来说，是早已解决了的问题：但由于其计算工作 量大，所以一直未能付诸实现。近年来，由于电子计算技术的飞速发展，为对电 场做快速和精确的数值模拟提供了技术条件，使之得以实现。目前，国内外已经 有不少单位和学者，成功地将数值模拟技术应用子研究激电异常。他们研究激电 异常的途径也是基于“等效电阻率法”原理i l 】，先模拟计算没有激电效应的一次 场，然后将地下各地质体的真电阻率换成相应的等效电阻率尼= 肛( 1 一叩，) ，并 6 中南大学硕士学位论文第2 章激发极化法及其数值模拟 计算二次场和极化率。所以实际上，这种对激电场的模拟方法还是归结为模拟不 包括激电效应场的一次场( 即稳定电流场) 。目前，用来对稳定电流场作数值模 拟的方法主要有两种，即有限差分法和有限单元法i s 。 2 2 1 有限差分法 有限差分法 3 9 1 是从电场( 电位) 或电磁场所满足的偏微分方程和边界条件 出发，将微分方程转变为差分方程，其研究步骤是：首先将研究区域按一定方式 离散化，然后在每个单元内设位、场呈线性变化，电性为均匀的，因而微分方程 的微分就可用差分来代替，于是就可以建立一组线性差分方程，最后求解此线性 方程组即得相应的位场分布。 。 通常解决稳定电流场空间分布的方法是根据电磁场理论导出一定的微分方 程，并用解析法求解；有限差分法则是将导出的微分方程离散化，得到相应的差 分方程以求解。具体来说，在稳定电流场中，电流密度_ ，、电场强度e 、电位v 和电导率口之间有如下关系： j = a e ，点= - v y ( 2 一i ) 即： ；：仃v 矿 ( 2 2 ) 在电流源以外的点上，根据电荷守恒定律有： v p - v 乃= 0 ( 2 3 ) 或者经过矢量运算，将上式改写为： v 2 ( 口矿) + ( ，v 2 矿一 ，v 2 盯= 0 ( 2 4 ) 上式便是在非均匀介质中，稳定电流场电位v 满足的微分方程。在均匀或分区域 均匀介质中，它简化为众所周知的拉普拉斯方程：v 2 = 0 。 在用有限差分法求微分方程( 2 - 4 ) 的数值解时，首先要在求解区( 或模拟 区) 建立网格，以便将微分方程离散化。对于模拟区( 即网格) 边缘上的节点( 外 边界点) ，其电位值可根据具体的地电条件给出，或者表示为该点坐标的已知函 数，或者表示为相邻节点电位的线性函数。这样，网格上的每一个节点都有一个 表示该点电位与相邻点电位之间关系得线性方程，总合起来，便得到一个未知数 ( 即各点的电位值) 和方程个数等于网格总节点数的高阶线性方程组。现代计算 技术有多种方法可以用来解这个方程组。由此，便可算出网格中各节点的电位值， 即求得模拟区内空间各点的电位分布。 7 中南大学硕士学位论文 第2 章激发极化法及其数值模拟 2 2 2 有限单元法 有限元法【舶】是从位、场所满足的偏微分方程出发，根据微分方程解与泛函 极小问题的等价性，将微分方程和其边界条件转化为相应泛函的变分阻题，其研 究的步骤是将研究区域按一定方式离散化，设单元内位、场呈线性变化，电性参 数均匀，这时泛函是各节点位、场的二次函数，利用求极小的必要条件，即泛函 对各节点位、场的变分为零，二次函数的交分为一次函数，由此得到个线性方 程组，解此线性方程组便可得各节点的位、场值。 有限单元法【9 】是近二十多年来随有限差分法之后迅速发展起来的又一种有效 的数值模拟方法。七十年代初，h j 科冈首先将其用于电法模拟。以后其他一些 学者进一步完善了它，使之成为求解复杂条件下电阻率法和激电法问题的有力工 具。 有限单元法在作激电法模拟时，也是分别求解地电断面为真电阻率和有效电 阻率分布时的稳定电流场( 即一次场和总场) ，然后根据等效电阻率法计算激电 异常。这与有限差分法的做法是一样的，所以，下面只讨论求解稳定电流场的电 位问题。 有限单元法求解稳定电流场电位v ，是根据变分法原理或最小总电能原理将 给定边界条件下求解( v ) 微分方程的问题，等价地转换成求v 的泛函的极值问 题；再经过离散化，得到由空间各点未知电位值组成的高阶线性方程组；最后在 数字电子计算机上求解该方程组，以确定空间各点的电位。 在包括电流源在内的不均匀三度空间中，稳定电流场的电位v 满足泊松方 程： v ( 口v n = 一v z ( 2 5 ) 式中盯为介质电导率；v 为稳定电流场电位：v 丘一电流源体密度。它们 都是空间坐标的函数。在直角坐标系中，公式( 2 5 ) 可写成： 丢p 篆，+ 茜p 号，+ 毫p 事= - v z c z 嘞 根据变分原理，若在区域q 的一部分边界黾上，v 已知；在其余边界屯上， v 满足“自然边界条件”：娶l ，其中n 是龟的外法线方向。则求解微分方程( 2 6 ) c 忡l 如 可等价转换为泛函极值问题，即找出一个函数v ，使其泛函j 取得极小值： 8 中南大学硕士学位论文第2 章激发极化法及其数值模拟 ，= 职竹警等2 + c 参2 书v z y 陋：= 以知；岛+ 是为q 的全部边界 = 极小 ( 2 7 ) 有限单元法和有限差分法是根据不同的原理作数值模拟，实践证明，两种方 法可以得到彼此相近的模拟结果。 有限差分法的程序编制比较简单，采用这种方法对平面规则形体( 如板状体、 层状体等) 进行数值模拟比较方便；但它对内节点、内边界点及两种类型的外边 界点要分别用不同的电位计算格式。所以，当不同电导率的介质较多或不均匀体 表面积为曲面时，使用这种方法不甚有利。 有限单元法易于模拟具有不规则表面的地电模型；此外，这种方法在整个求 解区可采用同一种计算格式，因而可避免采用多种计算格式时给模拟结果带来的 误差。但有限单元法程序编制较复杂，占用计算机的内存量也较大l l 】。 本文所用的地电模型主要是是有地形的模型，故采用有限元法进行实现。 a n s y s p l 是融结构、热、流体、电磁和声学为体的大型c a e 通用有限元软 件，可广泛地应用于核工业、铁道、石油化工、航空航天、机械制造、能源、汽 车交通、电子、土木工程、生物医学等一般工业及科学研究。该软件提供了不断 改进的功能清单，具体包括；结构高度非线性分析、电磁分析、计算流体力学分 析、设计优化、接触分析、自适应网络划分及利用a n s l 】s 参数设计语言扩展宏命 令等功能。另外，a n s y s 本身不仅具有较为完善的分析功能，同时为用户自己进 行二次开发提供了友好的开发环境，a p d l 语言和u i d l 语言是其开发的操作语言。 2 3a n s y s 软件介绍 a n s y s 的软件创始人是j o h ns w a n s o n 博士，在s w a n s o n 博士的领导下， a n s y s 软件的研究与开发不断汲取当今计算方法和计算机技术的最新发展，领 导着有限元界发展的趋势，并为全球工业界所广泛接受，拥有全球最大的用户群。 a n s y s 提供给用户的操作方式有图形用户界面g u i ( g r a p h i c a lu s e r i n t e r f a c e ) 、a n s y s 参数化设计编程语言 p d l ( a n s y sp a r a m e t e rd e s i g n l a n g u a g e ) 和用户界面设计语言u i d l ( u s e ri n t e r f a c ed e s i g nl a n g u a g e ) ，后 两者也是a n s y s 优化设计、自适应网格以及二次开发的主要基础。一般而言， a n s y s 的基本分析过程可以分为三步，即： 前处理( p r e d r o c e s s o r ) 施加荷载与求解( s o l u t i o n ) 后处理( p o s t p r o c e s s o r ) a n s y s 对每一个分析步骤都提供有定的工具模块。 9 中南大学硕士学位论文第2 章激发极化法及其数值模拟 前处理主要在前处理器( p r e p r o c e s s o r ，p r e p 7 ) 中完成，包括有限元模型 的创建等。前处理的主要任务是建立结构分析的有限单元模型，它是结构分析的 开始，一般步骤为嘲： ( 1 ) 分析准备 ( 2 ) 设置单元类型 ( 3 ) 设定实常数 ( 4 ) 定义材料属性 ( 5 ) 创建模型 ( 6 ) 划分网格 施加荷载与求解主要在求解器( s o l u t i o n ，s o l u ) 中完成，主要操作包括 定义自由度、施加荷载、设定分析类型以及求解等。施加荷载与求解的主要任务 是进行结构的计算，一般步骤为： ( 1 ) 施加荷载条件 ( 2 ) 给定边界条件及求解 后处理在后处理器中完成，a n s y s 有通用后处理器( g e n e r a l p o s t p r o c e s s o r ，p o s t l ) 和时程后处理器( t i m e h i s tp o s t p r o c e s s o r p s c r r 2 6 ) 两 种。前者主要用于静态分析的结果后处理，而后者主要用于动态分析结果的后处 理。后处理的主要任务是进行结构的分析。 应用a n s y s 进行电法有限元正演模拟分析。创建有限元模型是整个正演模 拟的基础，模型的建立直接关系到模拟的速度和精度，在建立模型时，我们可以 进行合理的简化，尽量使进行的网格划分和加载变得简单。在划分网格时，我们 可以在电源点和异常体附近划分较密的网格，随着与电源点的距离的增加，可以 把网格划的稀疏一点，这样在不增加总节点数的情况下，可以提高我们所关心区 域的网格划分的密度，从而提高精度、加快计算速度。下面着重介绍a n s y s 建 模与网格划分。 2 4a n s y s 有限元建模与网络划分 有限元分析是针对特定的模型而进行的，因此，用户必须建立一个有物理原 型的准确的数学模型。一般情况下，有限元模型包括所有的节点、单元、材料属 性、实常数、边界条件、荷载条件以及用这些条件来定义的物理场特征。建模与 网格划分是有限元分析的基础。 利用a n s y s 软件前处理模块的模型生成功能，可以快速、方便地建立实际工 程系统的有限元模型。 1 0 中南大学硕士学位论文第2 章激发极化法及其数值模拟 2 4 1 生成模型的典型步骤 通常的建模过程【4 1 1 如下： 一设计建模方案。首先确定分析目标，设计模型采取的基本形式，选择合 适的单元类型，并考虑如何建立适当的网格密度； 进入前处理器开始建立模型； 建立工作平面； 利用几何元素和布尔运算操作生成基本的几何形状； 激活适当的坐标系； n 生成单元属性表( 单元类型、实常数、材料属性和单元坐标系) ； 一设置属性指针； 通过对实体模型划分网格来生成节点和单元； 模型数据保存； 退出前处理。 2 4 2 d g s y s 有限元网格划分 现今所有的有限元分析都是采用实体建模，j 心i s y s 以数学的方式表达结构的 几何形状，用于在里面填充节点和单元，还可以在几何边界上方便地施加载荷， 但是几何实体模型并不参与有限元分析，所有施加在有限元边界上的载荷或约 束，最终必须传递到有限元模型上( 节点和单元) 进行求解。 通常网格划分有三个步骤 5 】： 建立选取单元数据 一设定网格建立所需的参数 一产生网格与修正网格 第一步是建立单元的数据，这些数据包括单元的种类( t 1 p e ) ，单元的几何 常数( r ) ，单元的材料性质( m p ) ，及单元形成时所在的坐标系，以及不同的单 元坐标系统。 第二步是设定网格划分的参数，最主要是定义对象边界单元的大小和数目。 网格设定所需的参数，将决定网格的大小和形状，这一步将影响分析结果的正确 性和过程的经济性。 完成前两步之后即可进行网格划分，即进行第三步的工作，完成有限元模型 的建立。如果不满意网格划分的结果，也可以清楚网格划分，重新定义单元的大 小和数目，再进行网格划分，直到得到满意的有限元结果为止。 有限元方法的基本思想是将结构离散化，用有限个容易分析的单元来表示复 中南大学硕士学位论文 第2 章激发极化法及其数值模拟 杂的对象，单元之间通过有限个节点相互连接，然后根据变形协调条件综合求解。 由此，有限元网格的划分一方面要考虑对各物体几何形状的准确描述，另一方面 也要考虑变形梯度的准确描述。网格划分的密度是个重要的问题，太密则会大大 增加计算时间，但计算精度却不会成比例地增加。这样，就存在一个最佳网格密 度问题，这个问题往往需要多年工作经验的积累。如果前置处理程序能够自动确 定网格密度，对节省机时的意义非常大。对于地球物理模型而言，一般都是一些 较复杂的地电模型，因此建立适当的地电模型和划分适当网格应遵循一定原则， 建立复杂有限元网格划分的基本原则如下冈： 1 ) 网格的大小 网格大小就是所取泛函积分区域的大小，一般来说网格越大越好。对于微分 方程边值问题的求解，只有给出正确的边界条件，才能求出域中比较精确的函数 值。例如，对于使用第一类边界条件时，要求给定正确的边界函数值。但是对于 不均匀的地电断面，边界函数值和其他边界条件值，特别是地下部分，均无法正 确求出。因此往往采用地下边界值为零，或由均匀地电条件给值。这时就要求区 域的边界远离不均匀区，即要求网格大，否则就会影响计算的精度，但是，另一 方面，如果网格内单元的大小不变的话，那么网格太大，势必要大量的增加节点 数，从而需要更多的计算机内存和增加计算工作量，这是因为在a n s y s 里面， 节点个数的增加，将使线性方程组的阶数增加，于是刚度系数矩阵的形成，解此 方程的时间必然会增加。对于网格内部单元的大小，一般来说单元越小，计算精 度越高。我们假定u 函数在每个单元内呈线性变化。如果单元太大，实际函数便 可能不满足这个条件，从而增加计算误差。我们还假定单元内电性是均匀的，即 电导率为常数，这也要求单元较小，特别是要拟合复杂的地电断面和地形剖面， 更需要划分得细致些，才能满足单元内电性均匀的条件。为了克服网格大小和单 元大小选择精度和工作量之间的矛盾，我们在划分网格时采用非均匀的网格，网 格的中心部分单元小，节点密，边界单元大，节点稀，由中心到边缘单元逐渐放 大，这样既保证了网格有足够的大小，又保证地电断面的复杂部位位于网格中心， 以满足单元内电性均匀和u 函数线性变化的条件。 2 ) 网格密度 有限元网格数量的多少将直接影响计算结果的精度和计算规模的大小。一般 来讲，网格数量增加，计算精度会有所提高，但同时计算规模也会增加，所以在 确定网格数量时应权衡这两个因素综合考虑。网格较少时增加网格数量可以使计 算精度明显提高，而计算时间不会有大的增加。当网格数量增加到一定程度后， 再继续增加网格时精度提高甚微，而计算时间划分的计算结果，如果两次计算结 果相差较大，应该继续增加网格，重新计算，直到误差在允许范围之内。在决定 1 2 中南大学硕士学位论文第2 章激发极化法及其数值模拟 网格密度时还应考虑分析类型。对于进行电磁场分析时( 对本文而言特指直流电 场分析) ，网格类型与源点的位置对计算精度结果是有影响的，在二维计算中， 网格形状的不对称，源点位置的不对称，都会引起计算结果的不对称：在地下半 空间中，源附近、近地表附近及异常体处的电位的变化比较剧烈，所以该处的网 格应较密，其它地方适当可以稀疏一些。 3 ) 单元阶次 a n s y s 软件中单元类型多走l t o o 多种，适合电磁场分析模块的单元类型就有 几十种，许多单元都具有线性、二次和三次等形式。其中二次和三次形式的单元 称为高阶单元。选用高阶单元可提高计算精度，因为高阶单元的曲线或曲面边界 能够更好地逼近结构地曲线和曲面边界，且高次插值函数可以更高精度地逼近复 杂场函数，所以当结构形状不规则、场的分布复杂时可以选用高阶单元。增加网 格数量和单元阶次都可以提高计算精度。因此在精度一定的情况下，用高阶单元 离散模型时应选择适当的网格数量，太多的网格并不能明显提高计算精度，反而 会使计算时间大大增加。为了兼顾计算精度和计算量，同一模型可以采用不同阶 次的单元，即精度要求高的重要部位用高阶单元，精度要求低的次要部位用低阶 单元 4 ) 网格质量 它是指网格几何形状的合理性。质量好坏将影响计算精度，质量太差的网格 甚至会中止计算。直观上看，网格各边或各个内角相差不大，网格面不过分扭曲， 边节点位于边界等份点附近的网格质量较好，网格质量可用细长比，锥度比、内 角、翘曲量、拉伸值、边节点位置偏差等指标度量，这些参数均可以利用a n s y s 单元检查功能获得。一般情况下，对于二维分析问题，单元形状应该取四边形； 在三维分析问题中，通常将单元分成规则六面体形状，这样可以有效地控制单元 和节点规模，并且能有效保证单元质量。 5 ) 网格剖分的自适应性 自适应网格划分是有限元分析中网格划分的主要思想方法之一。自适应是一 种自动近似的技术，程序可以自行减少由于网格划分不当引起的误差，而不需要 人为地指定单元的大小。自适应有限元分析的最基本思想就是通过事后误差估 计，对计算结果误差较大的区域进行局部网格优化。自适应有限元分析的网格划 分，要求不仅在生成初始网格时能完全自动有效，而且应具有极强的、面向整个 分析域的网格控制能力和局部网格修正能力。并且进行自适应有限元分析有一个 非常重要的前提条件，就是能在任意的几何分