ANSYS在直流电法正演中的应用

第21卷　第3期地　球　物　理　学　进　展Vol.21　No.32006年9月(页码:987～992)PRO GRESS IN GEOP H YSICSSept.　2006ANSYS 在直流电法正演中的应用汤井田,　肖　晓,　杜华坤,　王　武(中南大学信息物理工程学院,长沙410083)摘　要　以点电源场为例,详细叙述了直流电法正演在ANSYS 上的实现过程.通过正演模型ANSYS 模拟解与理论值的对比分析,验证了应用ANSYS 进行直流电法正演模拟的正确性.鉴于ANSYS 强大的分析功能,快速、准确的模拟以及便于操作等优点,对ANSYS 在地球物理正演中的应用的研究是有重大意义的.关键词　ANSYS ,正演,直流电法,有限元中图分类号　P631 文献标识码　A 文章编号　100422903(2006)022*******The Application of ANSYS in direct current method forw ard modelingTAN G Jing 2tian ,　XIAO Xiao ,　DU Hua 2kun ,　WAN G Wu(School in I nf o 2physics and Geomatics Engi neeri ng ,Cent ral S out h Universit y ,Changsha 410083,China )Abstract With point electric source as example ,The paper introduces process of carying out DC method forward modeling with ANSYS.Bycomparing and analyzing ANSYS modeling result and theoretic value of typical model ,It proves the validity of DC method forward modeing with ANSYS.ANSYS has the virtue for powerf ul analysis capabil 2ity ,quick and accurate forward modeing and convenient operation.It is important in appying ANSYS in geophysical modeling.K eyw ords ANSYS ,forwardmodeling ,DC method ,finite element method收稿日期　2005207210;　修回日期　2005208220.基金项目　国家自然科学基金重大项目(50099620203202)资助.作者简介　汤井田,男,汉,1965年生,博士,中南大学信息物理工程学院教授,长期从事应用地球物理研究工作,主要研究地球物理信号处理的现代数学方法、资源地球物理勘探、工程及灾害地球物理勘探和地球物理数据高分辨处理与解释等方面.已发表学术论文70余篇.(E 2mail :jttang @ )0　引　言ANS YS 软件是一个以有限元分析为基础的大型通用CA E 软件,它具有强大而广泛的分析功能.它包括热、电、磁、流体和结构等诸多模块,具有强大的求解器和前、后处理功能,为我们解决复杂、庞大的工程项目和致力于高水平的科研攻关提供了一个优良的工作环境,更使我们从繁琐、单调的常规有限元编程中解脱出来.ANS YS 本身不仅具有较为完善的分析功能,同时也为用户自己进行二次开发提供了友好的开发环境.ANS YS 软件自20世纪70年代以来,不断吸收新的计算方法和计算技术,其发展一直处于前列.对于直流电法,已知电阻率的空间分布求电场分布的过程称为正演.正演是我们进行地球物理反演和解决实际地球物理问题的基础,也是地球物理工作者比较关注的.本文将以点电源场为例,详细介绍应用ANSYS 的静电模块进行直流电法正演的过程.1　点电源场的基本方程在稳定电流场中,若电流密度为j ,电场强度为E ,电位和介质的电导率存在如下的关系[1～6]:j =σE 和E =-Δu .在地面A (x A ,y A ,z A )点有一电流强度为I 的点电源,用j 表示电流密度矢量,在空间作任意闭合面Γ,Ω是Γ所围的区域,则区域内的电位应满足微分方程[1,4]:Δ・(σΔu )=-2I δ(A ).(1)直流电法要解决的问题就是寻求方程(1)满足一定边界条件的解.然而只有在一些特殊情况下上式才有解析解,一般实际工作中常采用数值解法[7～21].地　球　物　理　学　进　展21卷图1　点电源边界示意图Fig.1　Sketch of point electricsource boundary 在直角坐标系中,点电源电位标量的边值问题可归纳为[1,4]:9 9x σ9u9x+99yσ9u9y+99zσ9u9z =-2Iδ(x A)δ(y A)δ(z A),(2)σ19u19n1+σ29u29n2Γ1=0,[u1-u2]Γ1=0,(3) 9u9nΓs=0,(4)u|Γ∞=cr.(5)式中ΓS表示地面边界,Γ1表示区域内的介质分界面,Γ∞表示无穷远边界(见图1).δ(x A)表为钬利特函数,σ1、σ2为内部分界两侧的电导率,r为Γ∞上点与A点的距离,v1v2为内部分界面后两侧电位.2　ANSYS模拟ANS YS提供给用户的操作方式有图形用户界面(Grap hical U ser Interface,GU I)和ANS YS参数化设计编程(ANS YS Parameter Design Lan2 guage,A PDL)两种形式,后者也叫命令流形式,两种形式互补,给用户带来了极大的方便.同时,后者也是ANS YS优化设计、自适应网格以及二次开发的主要基础.ANS YS分析可以分为三大步骤[22]:创建有限元模型、施加载荷并求解、后处理(见图2).创建有限元模型是整个正演模拟的基础,模型的建立直接关系到模拟的速度和精度.一般的地球物理模型相对来说都是比较简单的,所以有限元模型的创建也是很简单的.值得一提的是,在建立几何模型时,我们可以进行合理的简化和近似,尽量能使接下来将要进行的网格划分和加载变得简单.例如,我们在进行点电源半无限空间的三维模拟时,就可以考虑是不是可以用一个半球体来代替半无限空间,并把点电源置于球心,这样不管是网格划分还是加载都将变的容易控制.在划分网格时,考虑到ANS YS网格化的容易控制的特性,我们可以在电源点和异常体附近划分较密的网格,随着与电源点的距离的增加,我们把网格划的越来越稀疏.这样在不增加总节点数的情况下,可以提高我们所关心区域的网格划分的密度,从而提高精度、加快计算速度.图2　ANSYS分析步骤框图Fig.2　Chart of ANSYS analysis step 施加载荷与求解是进行ANS YS模拟的第二步,对于直流电法的正演模拟实际上就是确定供电电流和加边界条件的过程,这是整个模拟中最关键的一步.对于内边界和地表边界(3式和4式)在ANS YS计算过程中是自动满足的,可以不予考虑.对于无穷远边界,我们可以根据不同的装置进行加载.以偶极装置为例,每移动一次电极就需要进行重新加载和求解,我们可以用ANS YS提供的A PDL 语言进行此操作.命令行如下:3SET,Aa,2 !定义AB、MN的长度3SET,Mm,12 !定义计算范围(-8893期汤井田,等:ANS YS在直流电法正演中的应用Mm,Mm)3SET,Res1,500 !定义背景电阻率3SET,Res2,100 !定义异常区电阻率3SET,Currt,1 !定义电流强度3SET,Pi,ACOS(21) !定义参数Pi3DIM,V1,A RRA Y,(23Mm/Aa+1)35,,　!定义一维数值型数组3DIM,Ra,TABL E,(23Mm/Aa+1)35,2,　!定义二维表格型数组(存放电阻率及其对应坐标)3DO,Nn,1,5,1 !从n=1循环到n=53DO,I,1,23Mm/Aa+1,1 !每条剖面从2Mm到Mm以一米为间隔逐点计算/SOL U !进入求解器L SCL EA R,ALLNSEL,S,EXT !选择模拟无穷远边界上的节点NSEL,U,LOC,Y3GET,Nnod,NODE,,COUN T3GET,Nmin,NODE,,NUM,M IN3DO,J,1,Nnod,1 !逐点赋电位值Dist1=SQ R T((NX(Nmin)+(Nn/2+1)3Aa +Mm2(I21)3Aa)332+N Y(Nmin)332) Dist2=SQ R T((NX(Nmin)+(Nn/2+1)3Aa +Mm2I3Aa)332+N Y(Nmin)332)V0=Currt3Res13(1/Dist121/Dist2)/(23 Pi)D,Nmin,VOL T,V0Nmin=NDN EXT(Nmin)3ENDDONSEL,ALL !在A、B两点赋电流值NSEL,S,NODE,,NOD E(2(Nn/2+1)3Aa2 Mm+(I21)3Aa,0,0)F,ALL,AM PS,CurrtNSEL,ALLNSEL,S,NODE,,NOD E(2(Nn/2+1)3Aa2 Mm+I3Aa,0,0)F,ALL,AM PS,2CurrtALL SEL,ALLSOL V E !求解FIN I/POST1 !后处理,读取,并计算视电阻率V1(I+(Nn21)3(23Mm/Aa+1))=VOL T (NOD E(I3Aa2(12Nn/2)3Aa2Mm,0,0))2VOL T (NODE(Nn3Aa/22Mm+I3Aa,0,0))Ra(I+(Nn21)3(23Mm/Aa+1),1)=2Mm +(I21)3AaRa(I+(Nn21)3(23Mm/Aa+1),2)=V1(I +(Nn21)3(23Mm/Aa+1))3Pi/Currt/LO G ((Nn+2)3Nn/(Nn+1)332)3ENDDOFIN I3ENDDO/GCOL,1,Resa3V PLO T,,Ra(1,1),2,ANSYS的后处理功能是非常强大的,计算完毕我们可以很方便的查看节点或者单元的电位分布图,并可以读取任意节点或单元的电位并可对其进行简单的计算.如上面命令行中,我们很方便地由模拟出的电位计算出视电阻率.3　ANSYS模拟的误差分析为了验证ANS YS模拟结果的的正确性,我们对图3所示的二层地电断面进行了模拟值与理论值的对比.图中,X表示电源点A到M N中点的距离.我们首先用解析法计算了不同的X处的M、N之间的电位差(ΔU MN),如表1第二列所示;表1第三列为有限元模拟的结果.图3　二层地电断面示意图Fig.3　Skech of two2layers geolelectric cross section表1　点电源层地电断面电位差的有限元计算结果与解析解的对比Table.1　Comparing of finite element modeling resultand theoretic value of two2layers geoelctric crosssection with point eletric sourceX(m)理论值(v)模拟值(v)误差(v)1.522.261452822.24490410.01654872.513.233455913.23529060.00183473.59.609565619.616007190.006441584.57.675227177.684581660.00935455.56.48943311 6.50128900.01185598.5 4.73045518 4.749297450.018842310.5 4.17116734 4.194613260.023445912.5 3.80950693 3.837594460.0280875989地　球　物　理　学　进　展21卷图4　模型一的几何模型示意图Fig.4　Sketch of geometry model of model one 从表1看到,在电源点附近(小于1.5m )误差较大,从2.5m 以后,误差随着X 的增大而增大,但总体上来说绝对误差小于0.03伏,相对误差小于0.8%.这样的计算精度是令人满意的.4　计算实例模型一:水平大地下,有一长方体(长方体的长边沿Z 轴,长远远大于高与宽),其横截面为5m ×3m 的长方形,长方体的电阻率ρ2=100Ω・m ,埋深为2m ,大地电阻率ρ1=500Ω・m 根据模型的特征,可以进行二维模拟(在x y 平面内进行模拟),几何模型如图4.图5分别为中梯装置和三极装置的视电阻率曲线(在模型和网格划分都一样的情况下),它们都清楚地反应出一个位于-2.5m 到2.5m 之间的一个低阻异常.在小号点处两条曲线基本重合,但到大号点处三极装置的视电阻率曲线出现向上偏移的趋势,这是由于ANS YS 算出的电位随着距电源点距离的增加精度也有所降低所造成的.图5　正演视电阻率曲线Fig.5　Apparent restivity curve of forward modeling图6和图7分别是a =1和a =2时偶极装置视电阻率断面图,从图6可以看出位于中心位置的图6　偶极装置视电阻率断面图A B =M N =1,　N =1,2,3,4,5Fig.6　Apparent resistivitycross section with dipole 2dipole图7　偶极装置视电阻率断面图A B =M N =2,　N =1,2,3,4,5Fig.7　Apparent resistivity cross section with dipole 2dipole993期汤井田,等:ANS YS 在直流电法正演中的应用一低阻异常体顶部位于n =3左右,根据深度和n 经验转换公式[3]可知在对应深度约2m 左右,这与正演模型是吻合的.同时,图6不难看出在-2.5m 到2.5m 之间的异常顶部是一条直线,这也正好反映出长方形异常的宽度.从图7可以看出一个低阻异常位于n =1到n =4之间,对应的深度在2m 到5m 之间,这与正演模型吻合.图8　模型一的几何模型示意图Fig.8　Sketch of geometry m odel of m odel two 模型二:水平大地下,有一圆柱(圆柱的走向沿Z 轴),其横截面为r =2m 的圆,圆柱的电阻率ρ2=500Ω・m ,埋深为2m ,大地电阻率ρ1=100Ω・m.在一断面上进行二维模拟,建立如图8所示的几何模型.图9是在模型具有相同的网格划分时,进行的中梯和三极装置模拟的视电阻率曲线,两条曲线都很明显地反映出在-2m 到2ms 之间的高阻异常,其中,中梯装置的效果更好.图9　正演视电阻率曲线Fig.9　Apparent resistivity curve of forward modeling 图10和图11分别是a =1和a =2时偶极装置视电阻率断面图,从图10可以看出位于中心位置的一低阻异常体顶部位于n =3(深2m )左右,这与正演模型是吻合的.与图6相比,我们可知图10所反映的异常体顶部是弧形的,这与正演演模型中的圆形异常体是吻合的.图11显示,在点号为零处的电图10　偶极装置视电阻率断面图A B =M N =1,　N =1,2,3,4,5Fig.10　Apparent resistivity cross section with dipole 2dipole图11　偶极装置视电阻率断面图A B =M N =2,　N =1,2,3,4,5Fig.11　Apparent resistivity cross section with dipole 2dipole199地　球　物　理　学　进　展21卷阻率从随着n的增大从119Ω・m左右增大到134Ω・m,然后在减小到119Ω・m,相对于背景值100Ω・m来说,这是一个高阻区,它的深度范围是从2m到6m,这个异常范围和正演模型是相符的.5　结　论计算表明,用ANS YS进行直流电法的正演计算是可行的、可靠的.其强大的模拟功能和精确快速的计算能力为正演计算带来了许多便利.ANS YS 的A PDL语言可以灵活地实现有限元分析的众多功能,是进行ANS YS二次开发强有力的工具之一ANS YS的众多优越性表明它在地球物理的科研和教学中有着广阔的应用前景.参　考　文　献(References):[1]　徐世浙.地球物理中的有限单元法[M].北京:科学出版社,1994.[2]　屈超纯,张静,宋守根.点源场计算方法[M].昆明:云南科技出版社,1999.[3]　中南矿冶学院物探所教研室编.金属矿电法勘探[M].北京:冶金工业出版社,1980.[4]　黄俊革.三维电阻率/极化率有限元正演模拟与反演成像[D].[博士论文].长沙:中南大学,2001.[5]　周熙襄,钟本善,等.电法勘探数值模拟的若干结果[J].地球物理学报,1983,26(5):479～491.[6]　乔松,周锰钰,白朗.勘探电磁场论[M].北京:中国矿业大学出版社,1991.[7]　阮百尧,熊彬,徐世浙.三维地电断面电阻率测深有限元数值模拟[J]:地球科学,2001,26(1):73～77.[8]　阮百尧,徐世淅,等.三维地电断面电阻率测深有限元数值模拟[J].地球科学,2001,26(1):73～77.[9]　毛先进,鲍光淑.一种适于电阻率成像的正演新方法[J].地球物理学报,1998,41(增):385～393.[10]　底青云,倪大来,王若.高密度电阻率成像[J].地球物理学进展,2003,18(2):323～326.[11]　吴小平,汪彤彤.电阻率三维反演研究进展[J].地球物理学进展,2002,17(1):156～162.[12]　葛为中,层状介质点源电场正演解析及其应用[J].地球物理学报,1994,37(增):534～541.[13]　黄俊革,鲍光淑,阮百尧.坑道直流电阻率测深异常研究[J].地球物理学报,2005,48(1):222～228.[14]　阮百尧,熊彬.电导率连续变化的三维电阻率测深有限元模拟[J].地球物理学报,2002,45(1)131～138.[15]　Lamontagne Y,West G F.EM response of a rectangular t hinplate[J].Geophysics,1971,36:26.[16]　J epsen A F.Numerical modeling in resistivity prosperting[J].Ph. D.t hesis,University of california.Berkeley.1975.[17]　Mufti J R.Finite2difference resistivity modeling for arbitrari2ly shaped two2dimensional st ructures[J].Geophysics,1976,41:62.[18]　Mufti J R.A practical approach to finite2differ2ence modeling[J].Geophysics,1978,43:5.[19]　Dey A,Morrison H F.Resistivity modeling for arbitrarlyshaped t hree2dimensional struct ures[J].Geoppysics,1979,44:753.[20]　G ldman M M,Stoyert C H.Finite2difference calculations oft he transient field of an axially symmetric eart h for verticalmagnetic dipole excitation[J].Geophysics,1983,48:953. [21]　Spitzer K,Wurmstich B.Speed and Accuracy in32D resistici2ty modeling[J].Geophysics,2001,66.[22]　龚曙光,谢桂兰.ANSYS操作命令与参数化编程[M].北京:机械工业出版社,2004.299。