地球物理中的有限单元法

1、地球物理算法技术（论文）地球物理中的有限单元法院系:地球物理与信息技术院姓名:刘雅宁学号:201012005任课老师:张贵宾地球物理中的有限单元法一、有限单元法的介绍在地球物理理论计算中，存在着两类基本问题：正问题和反问题。给定场源的分布，求解场值的大小，这是正问题，或者称为正演问题。地球物理正演的数值计算方法，种类很多，最常用的有：有限差分法和有限单元法。有限单元法是50年代首先在弹性力学中发展起来的方法。主要优点是，适用于物性参数复杂分布的区域，但计算量大。随着计算机技术的发展，有限单元法在解决各个工程领域的许多数学物理问题中，得到了广泛的应用，称为一种高效、通用的计算方法。地球物理中的一

2、些边值问题，也采用了有限单元法，解决了许多从前无法计算的地球物理问题。有限单元法解决数学物理边值问题的基本思路和过程如下：1、给出地球物理边值问题中的偏微分方程和边界条件（及初始条件）。这一点看起来似乎容易，但做起来并不容易，特别是边界条件的给定。只有对地球物理方法的原理和问题有深入的理解，才能给边值问题中的偏微分方程和边界条件以正确的描述。2、将地球物理边值问题转变为有限元方程。实现这种转变的主要数学工具是变分法，用变分法得到的有限元法方程称为泛函极值问题。3、用优先单元法解决泛函极值问题其步骤大致如下：把研究区域剖分成有限个小单元，在每个单元上，把函数简化成线性函数、二次函数或高次函数，这

3、称为单元上函数的插值。用简化后的函数计算每个单元上的泛函。各单元之间，通过单元间节点上的函数值相互联系起来。对各单元的泛函求和，获得整个区域上的泛函。这样，有限单元法将连续函数的泛函，离散成各单元节点上函数值得泛函。根据泛函取极值的条件，得到各节点的函数值应满足的线性代数方程组。解代数方程组，得到各节点的函数值。有限单元法的主要优点是，适用于物性复杂分布的地球物理问题，而且，其解题过程也比较规范化。这些优点是有限单元法在地球物理中获得广泛的应用。但是，有限单元法是区域性方法，必须在全区域中进行剖分。剖分后的单元和节点数目多，最后得到的线性代数方程组很大。特别是三维问题和地球物理中常遇到的无解区

4、域问题，需要中、大型计算机，才能完成有限单元法的计算。这是有限单元法的主要缺点。e基本步骤用三角单元对区域进行剖分：用三角单元对整个区域进行剖分，因为三角单元的边容易拟合地形线的形状。三角形的顶点称为节点，用节点上的离散场值来近似场值的连续分布。剖分后对单元和节点进行编号，次序号（i,j,m）按逆时针方向排列。第一类边界条件的节点上的场值是已知的，其余节点上的场值是待求的；线性插值：在三角单元内假定场值u是线性变化的，则u=ax+by+c，u还可表示为u=Nu+Nu+Nu，其中NNN是形函数，iijjmm，ijmN一(ax+by+c)，N一(ax+by+c)，i2iiij2jjj1N一(ax+

5、by+c)，m2mmma=y-y,b=x-x,c=xy-xy，ijmimjijmmja=y-y,b=x-x,c=xy-xy，jmijimjmiima=y-y,b=x-x,c=xy-xy，mijmjimijji1一(ab-ab）是三角兀面积2ijji单元分析：将全区域的积分分解为单元e的积分之和F(u)=J(Vu)2dQ=2丄(Vu)dQ=21e翌)22,x1e+(聖)2dxdy,y继续推得单元e上的积分1dxdyquTku2eee其中uT”（uuu）是单元的u值列向量；k是单元系数矩阵,eijmek=丄（aa+b）bS,t=i，j，m，k是对称矩阵。st4stste总体合成:将单元的场值列向量

6、ue扩展成全体节点的场值向量u，口=（片口2.吟.Uj.um.und）,按照节点的总体序号，将单元系数矩阵中的各元放在k的相应行与列的交叉位置上，其余位置的元为零，这样单元积分可写成1Fe(U)2UTkeUe1uTku2e各单元积分相加时，只要将ke相加即可；求变分：通过以上四个步骤，已经将连续函数g的泛函，离散成各节点g值的多元函数：(u=-uTku，泛函的极值等于多元函数的极值，用多元函数求极值2的方法，对上式求微分，11,F(u)d(uT)=(duTku+uTkdu)=22因为K是对称矩阵，有duT=uTkdu，所以,F(u)duTku=0，由于duT丰0，所以由上式得ku=0。得到含有

7、ND个元的ND个方程联立的线性代数方程组；解线性代数方程组：首先将第一类边界条件代入，通过定带宽储存的对称带型线性方程组，解方程组，得到各节点的u值，至此有限单元法的求解过程结束。三、实现过程为了验证有限单元法的有效性，我们设计一个规则形状的地下矿体，给出模型：1、模型密度均匀的水平半无限空间，一个均匀球体S,球体半径R=10m,剩余密度a=1g/cm3，球心坐标(a，b)=(200,-100)。对于均匀球体来说，它与将其全部剩余质量集中在球心处的点的质量所引起的异常完全一样。若球心的埋藏深度为D，球的半径为R，剩余密度为a，则它的剩余质量为M=(4nR3a)/3,通过原点的任意水平剖面上则重

8、力异常的解析解表达式为：GMDg=(x2D2)3/2设测线长400m，高程变化(-200,300)，地形设为一曲线：y=20*sin(0.02*x)+30，其中x为测线离原点的水平距离，y为高程，则S引起的重力异常为：4冗R3gG(yb)3(xa)2(yb)2)3/2g=一2、剖分通过Matlab建模，得到地形曲线图，再用三角单元对划出的区域进行剖分(程序附后)，剖分后对单元和节点进行编号，并将节点的xy坐标和单元节点号列表(表1和表3),分别放在XY(2，ND)和13(3,NE)两个二维数组中。剖分图如下：内部的节点值。编制计算程100图吐：剖分图.0根据重力异常的解析解表达式算出区域边界上

9、及区域序，根据边界上节点的场值，算出区域内部节点的场值，将计算出来的内部节点场值与解析解比较，在有效数字内，若计算值与理论值一致或相差不大，则验证了有限单元法的可效性。3、剖分后的节点分布及单元编号见下表：节点总数ND=33；单元总数NE=40；3)单元节点编号表1单元序号12345678910111213141516171819i555588881111111114141414171717j123645697891210111215131415m41237456107891310111216131420212223242526272829303132333435363738394017202

10、02020232323232626262629292929323232321816171821192021242223242725262730282930331519161718221920212522232428252627312829304)节点坐标表2节点x000404040808080120120y3016530044.3172.130049.9174.930043.5171.747127355914795730927546412016016016020020020024024024028030028.83252164.4162630014.86395157.4319830010.07

11、671155.0383630017.37467280280320320320360360360400400400158.6873330032.33098166.1655030045.87336172.9366830049.78717174.89359300(5)第一类边界节点数ND1=24；6)第一类边界节点号和场值：表3边界节点序号123467910121315边界节2.6762.0231.2494.0311.3985.9141.5349.0421.64614.661.720点场值(*103)82048112854052080970478691587687926769839846916181

12、92122242527283031323321.181.74619.141.72011.441.6466.4871.5344.0161.3982.6831.9551.2492479724094638469563392676173915854140970269151458540根据书中给出的第一类边界条件、三角元剖分、线性插值的位场延拓得有限单元程序框图：编制计算程序，并将已知的边界节点数据输入，来计算区域内部的节点值（5、8、11、14、17、20、23、26、29点）。运行程序后，得计算结果如下：内部节点号数值解解析解误差50.00283778370.00241706500.00042071

13、8780.00364900470.00284539490.0008036098110.00450467130.00334078820.0011638831140.00540201320.00386391880.0015380944170.00614080070.00421715250.0019236482200.00609677100.00414288320.0019538878230.00496956430.00364161000.0013279542260.00371221450.00298882000.0007233946290.00273859080.00240874780.00032

14、98430小结通过计算出的解与用解析公式得到的解进行比较，误差相差不大，证明了有限元方法是一种有效的正演方法。附录：计算机程序地形及图形剖分程序：clc,clearx=0:400;y=-200:0;XY=meshgrid(x,y);y1=20*sin(0.02*x)+30;%地形%Model=50*exp(-(200-X).人2/120-(-100-Y).人2/120);%模型x1=40*x(1:11);y2=20*sin(0.02*x1)+30;%观测点%绘图figureplot(x1,y2,Rv,MarkerFaceColor,r,MarkerSize,8)legend(剖分图)holdo

15、narea(x,y1,linestyle,none)holdoncontourf(X,Y,Model,400,linestyle,none)set(gca,ticklength,0.00.0,FontName,Verdana,FontWeight,Bold,fontsize,12)%colormapjet(32)forj=1:length(x1)fori=1:3X1(i,j)=x1(j);endY1(1,j)=300;Y1(2,j)=300/2+y2(j)/3;Y1(3,j)=y2(j);endholdonplot(40*x(1:11),y2,Rv,MarkerFaceColor,r,Mark

16、erSize,8)legend(剖分图)fori=1:3holdonplot(x1,Y1(i,:),ko,x1,Y1(i,:),k,linewidth,1.2);endforj=1:length(x1)holdonplot(zeros(3,1)+x1(j),Y1(:,j),ko,zeros(3,1)+x1(j),Y1(:,j),k,linewidth,1.2);endfori=1:length(x1)-2holdonplot(x1(i:i+1),Y1(1,i),Y1(2,i+1),k,linewidth,1.2);holdonplot(x1(i:i+1),Y1(3,i),Y1(2,i+1),k

17、,linewidth,1.2);endholdonplot(x1(length(x1)-1:length(x1),Y1(1,length(x1)-1),Y1(2,length(x1),k,linewidth,1.2);holdonplot(x1(length(x1)-1:length(x1),Y1(3,length(x1)-1),Y1(2,length(x1),k,linewidth,1.2);k=0;forj=1:length(x1)fori=3:-1:1k=k+1;text(X1(i,j)+2,Y1(i,j)+10,num2str(k),color,r,fontsize,10,fontna

18、me,华文中宋);endendaxis(0400-200400)节点坐标值计算及比较的Fortran程序：PROGRAMsecondDIMENSIONX(1:3),Y(1:3),B(1:3),C(1:3)DIMENSIONXX(1:3,1:11),YY(1:3,1:11),XY(1:2,1:33)DIMENSIONDG(1:3,1:11)DIMENSIONUO(40),NDB(9)INTEGER:ND,NE,ND1,IEREAL,ALLOCATABLE:KE(:,:),SK(:,:),A(:,:),I3(:,:)REAL,ALLOCATABLE:P(:),NB1(:),U1(:),U(:)ND

19、=33NE=40ND1=24ALLOCATE(U(1:ND),P(1:ND),NB1(1:ND1),U1(1:ND1),I3(1:3,1:NE)!地质体的模型参数PI=3.14159!PI常量G=6.672E-11!万有引力常量Xl=200.;Yl=-100.!球心S坐标(xl,y2)R=10.!半径Pl=l.!密度DOI=l,llXX(1,I)=40.*(I-1)!观测点坐标xxYY(1,I)=20.*sin(0.02*XX(1,I)+30.!观测点坐标yyXX(2,I)=XX(1,I)!网格点坐标xx1YY(2,I)=(300.+YY(1,I)/2.!网格点坐标yy1XX(3,I)=XX(

20、1,I)!网格点坐标xx2YY(3,I)=300.!网格点坐标yy2ENDDO!计算异常场值DOI=1,3DOJ=1,11!异常体S异常值DG(I,J)=4.*PI*R*3.*P1*G*(YY(I,J)-Y1)&/(3.*(XX(I,J)-X1)*2.+(YY(I,J)-Y1)*2.)*(3.0/2.0)*1.0E9WRITE(*,*)DG(I,j)ENDDOENDDO!INPUTDATADOI=1,11K=3*(I-1)+1XY(1,K)=XX(1,I)XY(1,K+1)=XX(2,I)XY(1,K+2)=XX(3,I)XY(2,K)=YY(1,I)XY(2,K+1)=YY(2,I)XY(2

21、,K+2)=YY(3,I)U(K)=DG(1,I)UO(K)=DG(1,I)U(K+1)=DG(2,I)UO(K+1)=DG(2,I)U(K+2)=DG(3,I)UO(K+2)=DG(3,I)ENDDO!输入参数I3,存放单元节点编号OPEN(2,FILE=I3.txt)READ(2,*)DOI=1,NEREAD(2,*)I3(1,I),I3(2,I),I3(3,I)ENDDOCLOSE(2)WRITE(*,*)SUCCESS:InputI3.txtJ=0NDB=(/5,8,11,14,17,20,23,26,29/)DO111I=1,NDDOK=1,9IF(I.EQ.NDB(K)GOTO11

22、1ENDDOJ=J+1NB1(J)=IU1(J)=U(I)111CONTINUE!CALLFUNCTION函数!CallMBWCALLMBW(NE,I3,IW)WRITE(*,*)SUCCESS:CallMBWALLOCATE(KE(1:3,1:3),SK(1:ND,1:IW),A(1:ND,1:IW)!CallUK1CALLUK1(ND,NE,IW,I3,XY,SK)WRITE(*,*)Success:CallUK1!CallUB1CALLUB1(ND1,NB1,U1,ND,IW,SK,U)WRITE(*,*)Success:CallUB1!CallLDLTCALLLDLT(SK,ND,IW

23、,U,IE)WRITE(*,*)Success:CallLDLT!NB1U1OPEN(3,FILE=NB1U1.txt)WRITE(3,(2A15)NB1,U1DOI=1,ND1WRITE(3,(I15,f15.5)NB1(I),U1(I)ENDDOCLOSE(3)WRITE(*,*)Success:OutputNB1U1.txt!XYOPEN(3,FILE=XY.txt)WRITE(3,(2A15)X,YDOI=1,NDWRITE(3,(2F15.5)XY(1,I),XY(2,I)ENDDOCLOSE(3)WRITE(*,*)SUCCESS:OutputXY.txt!ResultOPEN(4

24、,FILE=Result.txt)WRITE(4,(4A15)节点号,数值解,解析解,误差DOI=1,NDWRITE(4,(I15,f15.10,f15.10,f15.10)I,U(I),&UO(I),U(I)-UO(I)ENDDOCLOSE(4)WRITE(*,*)Success:OutputResult.txt!DEALLOCATEDEALLOCATE(I3,KE,A,U,P,NB1,U1)ENDPROGRAMsecond!SUBROUTINE!MBW计算总体系数矩阵的半带宽SUBROUTINEMBW(NE,I3,IW)INTEGERMREALI3(1:3,1:NE)IW=0DOI=1,N

25、EM=MAX(ABS(I3(1,I)-I3(2,I),ABS(I3(2,I)-I3(3,I),&ABS(I3(3,I)-I3(1,I)IF(M+1).GT.IW)IW=M+1ENDDOEND!UK1集成总体矩阵SUBROUTINEUK1(ND,NE,IW,I3,XY,SK)REALI3(1:3,1:NE),XY(1:2,1:ND),SK(1:ND,1:IW)REALX(1:3),Y(1:3)REALKE(1:3,1:3)DOI=1,NDDOJ=1,IWSK(I,J)=0ENDDOENDDODOL=1,NEDOJ=1,3I=I3(J,L)X(J)=XY(1,I)Y(J)=XY(2,I)ENDDOCALLUKE1(X,Y,KE)DO1J=1,3NJ=I3(J,L)DO1K=1,JNK=I3(K,L)IF(NJ.LT.NK)GOTO11NK=(NK-NJ+IW)SK(NJ,NK)=SK(NJ,NK)+KE(J,K)GOTO111NJ=(NJ-NK+IW)SK(NK,NJ)=SK(NK,NJ)+KE(J,K)NJ=NJ+NK-IW1CONTINUEENDDORETURNEND!UKE1计算单元系数矩阵KESUBROUTINEUKE1(X,Y,KE)REALX(1:3),Y(1:3),C