地球物理中的有限单元法解剖.pdf

地球物理算法技术（论文）地球物理中的有限单元法 院 系:地球物理与信息技术院 姓 名:刘雅宁 学 号:2010120053 任课老师:张贵宾 地球物理中的有限单元法 一、有限单元法的介绍 在地球物理理论计算中，存在着两类基本问题：正问题和反问题。给定场源的分布，求解场值的大小，这是正问题，或者称为正演问题。地球物理正演的数值计算方法，种类很多，最常用的有：有限差分法和有限单元法。有限单元法是 50 年代首先在弹性力学中发展起来的方法。主要优点是，适用于物性参数复杂分布的区域，但计算量大。随着计算机技术的发展，有限单元法在解决各个工程领域的许多数学物理问题中，得到了广泛的应用，称为一种高效、通用的计算方法。地球物理中的一些边值问题，也采用了有限单元法，解决了许多从前无法计算的地球物理问题。有限单元法解决数学物理边值问题的基本思路和过程如下：1、给出地球物理边值问题中的偏微分方程和边界条件（及初始条件）。这一点看起来似乎容易，但做起来并不容易，特别是边界条件的给定。只有对地球物理方法的原理和问题有深入的理解，才能给边值问题中的偏微分方程和边界条件以正确的描述。2、将地球物理边值问题转变为有限元方程。实现这种转变的主要数学工具是变分法，用变分法得到的有限元法方程称为泛函极值问题。3、用优先单元法解决泛函极值问题其步骤大致如下：把研究区域剖分成有限个小单元，在每个单元上，把函数简化成线性函数、二次函数或高次函数，这称为单元上函数的插值。用简化后的函数计算每个单元上的泛函。各单元之间，通过单元间节点上的函数值相互联系起来。对各单元的泛函求和，获得整个区域上的泛函。这样，有限单元法将连续函数的泛函，离散成各单元节点上函数值得泛函。根据泛函取极值的条件，得到各节点的函数值应满足的线性代数方程组。解代数方程组，得到各节点的函数值。有限单元法的主要优点是，适用于物性复杂分布的地球物理问题，而且，其解题过程也比较规范化。这些优点是有限单元法在地球物理中获得广泛的应用。但是，有限单元法是区域性方法，必须在全区域中进行剖分。剖分后的单元和节点数目多，最后得到的线性代数方程组很大。特别是三维问题和地球物理中常遇到的无解区域问题，需要中、大型计算机，才能完成有限单元法的计算。这是有限单元法的主要缺点。二、基本步骤 用三角单元对区域进行剖分：用三角单元对整个区域进行剖分，因为三角单元的边容易拟合地形线的形状。三角形的顶点称为节点，用节点上的离散场值来近似场值的连续分布。剖分后对单元和节点进行编号，次序号（i，j，m）按逆时针方向排列。第一类边界条件的节点上的场值是已知的，其余节点上的场值是待求的；线性插值：在三角单元内假定场值 u 是线性变化的，则 u=ax+by+c，u 还可表示为iijjmmu=u+u+uNNN，其中ijmN N N是形函数，iiiijjjjmmmmijmimjijmmjjmijimjmiimmijmjimijjiijji11=a x+b y+c=a x+b y+c221=a x+b y+c2a=y-yb=x-xc=x y-x ya=y-yb=x-xc=x y-x ya=y-yb=x-xc=x y-x y1=a b-a b2NNN（），（），（），（）是三角元面积 单元分析：将全区域的积分分解为单元 e 的积分之和 222211ee111gg()=d=d=()+()dxdy222xyNENEF uuu（）（），继续推得单元 e 上的积分 22e11()=22TeeeeuuF udxdyu k uxy 其 中eijm()Tuuu u是 单元 的 u 值 列向量；ek是单元系 数矩阵，ststst1k=a a+b b4（），s，t=i，j，m，ek是对称矩阵。总体合成:将单元的场值列向量 ue扩展成全体节点的场值向量 u，u=(u1u2uiujumund)，按照节点的总体序号，将单元系数矩阵中的各元放在ek的相应行与列的交叉位置上，其余位置的元为零，这样单元积分可写成 11()22TTeeeeeF uu k uu k u 各单元积分相加时，只要将 ke相加即可；求变分：通过以上四个步骤，已经将连续函数 g 的泛函，离散成各节点 g 值的多元函数：1=u ku2TF（u），泛函的极值等于多元函数的极值，用多元函数求极值的方法，对上式求微分，11()d=du ku+u kdu=22TTTF uu（ku）（）0 因为 K 是对称矩阵，有du=u kduTTku，所以()du ku=0TF u，由于du0T，所以由上式得ku=0。得到含有 ND 个元的 ND 个方程联立的线性代数方程组；解线性代数方程组：首先将第一类边界条件代入，通过定带宽储存的对称带型线性方程组，解方程组，得到各节点的u 值，至此有限单元法的求解过程结束。三、实现过程 为了验证有限单元法的有效性，我们设计一个规则形状的地下矿体，给出模型：1、模型 密度均匀的水平半无限空间，一个均匀球体 S，球体半径 R=10m，剩余密度=1g/cm3，球心坐标（a，b）=（200，-100）。对于均匀球体来说，它与将其全部剩余质量集中在球心处的点的质量所引起的异常完全一样。若球心的埋藏深度为 D，球的半径为 R，剩余密度为，则它的剩余质量为 M=(4R3)/3,通过原点的任意水平剖面上则重力异常的解析解表达式为：2/322)(DxGMDg 设 测 线 长400m，高 程 变 化（-200,300），地 形 设 为 一 曲 线：y=20*sin(0.02*x)+30，其中 x 为测线离原点的水平距离，y 为高程，则 S 引起的重力异常为：3223/24(b)3(a)(b)RG ygxy 2、剖分 通过 Matlab 建模，得到地形曲线图，再用三角单元对划出的区域进行剖分（程序附后），剖分后对单元和节点进行编号，并将节点的 xy 坐标和单元节点号列表（表 1 和表 3），分别放在 XY（2，ND）和 I3（3，NE）两个二维数组中。剖分图如下：图 1：剖分图 根据重力异常的解析解表达式算出区域边界上及区域内部的节点值。编制计算程序，根据边界上节点的场值，算出区域内部节点的场值，将计算出来的内部节点场值与解析解比较，在有效数字内，若计算值与理论值一致或相差不大，则验证了有限单元法的可效性。3、剖分后的节点分布及单元编号见下表：（1）节点总数 ND=33；（2）单元总数 NE=40；（3）单元节点编号 表 1 单元 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 i 5 5 5 5 8 8 8 8 11 11 11 11 14 14 14 14 17 17 17 j 1 2 3 6 4 5 6 9 7 8 9 12 10 11 12 15 13 14 15 m 4 1 2 3 7 4 5 6 10 7 8 9 13 10 11 12 16 13 14 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 17 20 20 20 20 23 23 23 23 26 26 26 26 29 29 29 29 32 32 32 32 18 16 17 18 21 19 20 21 24 22 23 24 27 25 26 27 30 28 29 30 33 15 19 16 17 18 22 19 20 21 25 22 23 24 28 25 26 27 31 28 29 30（4）节点坐标 表 2 节点 1 2 3 4 5 6 7 8 9 10 11 x 0 0 0 40 40 40 80 80 80 120 120 y 30 165 300 44.3 4712 172.1 7355 300 49.9 9147 174.9 9573 300 43.5 0927 171.7 5464 12 13 14 15 16 17 18 19 20 21 22 120 160 160 160 200 200 200 240 240 240 280 300 28.8 3252 164.4 1626 300 14.8 6395 157.4 3198 300 10.0 7671 155.03836 300 17.3 7467 23 24 25 26 27 28 29 30 31 32 33 280 280 320 320 320 360 360 360 400 400 400 158.6 8733 300 32.3 3098 166.1 6550 300 45.8 7336 172.9 3668 300 49.7 8717 174.8 9359 300 （5）第一类边界节点数 ND1=24；（6）第一类边界节点号和场值：表 3 边 界 节点序号 1 2 3 4 6 7 9 10 12 13 15 边 界 节点场值（*103）2.676 8204 2.023 8112 1.249 8540 4.031 5208 1.398 0970 5.914 4786 1.534 9158 9.042 7687 1.646 9267 14.66 69839 1.720 8469 16 18 19 21 22 24 25 27 28 30 31 32 33 21.182479 1.7467240 19.149463 1.7208469 11.445633 1.6469267 6.4876173 1.5349158 4.0165414 1.3980970 2.6832691 1.9555145 1.2498540 根据书中给出的第一类边界条件、三角元剖分、线性插值的位场延拓得有限单元程序框图：编制计算程序，并将已知的边界节点数据输入，来计算区域内部的节点值（5、8、11、14、17、20、23、26、29 点）。运行程序后，得计算结果如下：内部节点号 数值解 解析解 误差 5 0.0028377837 0.0024170650 0.0004207187 8 0.0036490047 0.0028453949 0.0008036098 11 0.0045046713 0.0033407882 0.0011638831 14 0.0054020132 0.0038639188 0.0015380944 17 0.0061408007 0.0042171525 0.0019236482 20 0.0060967710 0.0041428832 0.0019538878 23 0.0049695643 0.0036416100 0.0013279542 26 0.0037122145 0.0029888200 0.0007233946 29 0.0027385908 0.0024087478 0.0003298430 小结 通过计算出的解与用解析公式得到的解进行比较，误差相差不大，证明了有限元方法是一种有效的正演方法。输入原始数据：节点总数 ND,单元总数 NE,第一类边界点数 ND1,节点坐标 XY(2,ND),单元节点编号 I3(3,NE),第一类边界节点号 NB1(ND1),第一类边界节点场值 U1（ND1）计算半带宽 形成总体系数矩阵 代入第一类边界条件 解线性方程组 输出数值解 附录：计算机程序 地形及图形剖分程序：clc,clear x=0:400;y=-200:0;X Y=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;%观测点%绘图 figure plot(x1,y2,Rv,MarkerFaceColor,r,MarkerSize,8)legend(剖分图)hold on area(x,y1,linestyle,none)hold on contourf(X,Y,Model,400,linestyle,none)set(gca,ticklength,0.0 0.0,FontName,Verdana,FontWeight,Bold,fontsize,12)%colormap jet(32)for j=1:length(x1)for i=1:3 X1(i,j)=x1(j);end Y1(1,j)=300;Y1(2,j)=300/2+y2(j)/3;Y1(3,j)=y2(j);end hold on plot(40*x(1:11),y2,Rv,MarkerFaceColor,r,MarkerSize,8)legend(剖分图)for i=1:3 hold on plot(x1,Y1(i,:),ko,x1,Y1(i,:),k,linewidth,1.2);end for j=1:length(x1)hold on plot(zeros(3,1)+x1(j),Y1(:,j),ko,zeros(3,1)+x1(j),Y1(:,j),k,linewidth,1.2);end for i=1:length(x1)-2 hold on plot(x1(i:i+1),Y1(1,i),Y1(2,i+1),k,linewidth,1.2);hold on plot(x1(i:i+1),Y1(3,i),Y1(2,i+1),k,linewidth,1.2);end hold on plot(x1(length(x1)-1:length(x1),Y1(1,length(x1)-1),Y1(2,length(x1),k,linewidth,1.2);hold on plot(x1(length(x1)-1:length(x1),Y1(3,length(x1)-1),Y1(2,length(x1),k,linewidth,1.2);k=0;for j=1:length(x1)for i=3:-1:1 k=k+1;text(X1(i,j)+2,Y1(i,j)+10,num2str(k),color,r,fontsize,10,fontname,华文中宋);end end axis(0 400-200 400)节点坐标值计算及比较的 Fortran 程序：PROGRAM second DIMENSION X(1:3),Y(1:3),B(1:3),C(1:3)DIMENSION XX(1:3,1:11),YY(1:3,1:11),XY(1:2,1:33)DIMENSION DG(1:3,1:11)DIMENSION UO(40),NDB(9)INTEGER:ND,NE,ND1,IE REAL,ALLOCATABLE:KE(:,:),SK(:,:),A(:,:),I3(:,:)REAL,ALLOCATABLE:P(:),NB1(:),U1(:),U(:)ND=33 NE=40 ND1=24 ALLOCATE(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 !万有引力常量 X1=200.;Y1=-100.!球心 S 坐标(x1,y2)R=10.!半径 P1=1.!密度 DO I=1,11 XX(1,I)=40.*(I-1)!观测点坐标 xx YY(1,I)=20.*sin(0.02*XX(1,I)+30.!观测点坐标 yy XX(2,I)=XX(1,I)!网格点坐标 xx1 YY(2,I)=(300.+YY(1,I)/2.!网格点坐标 yy1 XX(3,I)=XX(1,I)!网格点坐标 xx2 YY(3,I)=300.!网格点坐标 yy2 ENDDO !计算异常场值 DO I=1,3 DO J=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.0E9 WRITE(*,*)DG(I,j)ENDDO ENDDO !INPUT DATA DO I=1,11 K=3*(I-1)+1 XY(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,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,*)DO I=1,NE READ(2,*)I3(1,I),I3(2,I),I3(3,I)ENDDO CLOSE(2)WRITE(*,*)SUCCESS:Input I3.txt J=0 NDB=(/5,8,11,14,17,20,23,26,29/)DO 111 I=1,ND DO K=1,9 IF(I.EQ.NDB(K)GOTO 111 ENDDO J=J+1 NB1(J)=I U1(J)=U(I)111 CONTINUE !CALL FUNCTION 函数 !Call MBW CALL MBW(NE,I3,IW)WRITE(*,*)SUCCESS:Call MBW ALLOCATE(KE(1:3,1:3),SK(1:ND,1:IW),A(1:ND,1:IW)!Call UK1 CALL UK1(ND,NE,IW,I3,XY,SK)WRITE(*,*)Success:Call UK1!Call UB1 CALL UB1(ND1,NB1,U1,ND,IW,SK,U)WRITE(*,*)Success:Call UB1!Call LDLT CALL LDLT(SK,ND,IW,U,IE)WRITE(*,*)Success:Call LDLT!NB1U1 OPEN(3,FILE=NB1U1.txt)WRITE(3,(2A15)NB1,U1 DO I=1,ND1 WRITE(3,(I15,f15.5)NB1(I),U1(I)ENDDO CLOSE(3)WRITE(*,*)Success:Output NB1U1.txt!XY OPEN(3,FILE=XY.txt)WRITE(3,(2A15)X,Y DO I=1,ND WRITE(3,(2F15.5)XY(1,I),XY(2,I)ENDDO CLOSE(3)WRITE(*,*)SUCCESS:Output XY.txt !Result OPEN(4,FILE=Result.txt)WRITE(4,(4A15)节点号,数值解,解析解,误差 DO I=1,ND WRITE(4,(I15,f15.10,f15.10,f15.10)I,U(I),&UO(I),U(I)-UO(I)ENDDO CLOSE(4)WRITE(*,*)Success:Output Result.txt!DEALLOCATE DEALLOCATE(I3,KE,A,U,P,NB1,U1)END PROGRAM second !SUBROUTINE !MBW 计算总体系数矩阵的半带宽 SUBROUTINE MBW(NE,I3,IW)INTEGER M REAL I3(1:3,1:NE)IW=0 DO I=1,NE M=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+1 ENDDO END !UK1 集成总体矩阵 SUBROUTINE UK1(ND,NE,IW,I3,XY,SK)REAL I3(1:3,1:NE),XY(1:2,1:ND),SK(1:ND,1:IW)REAL X(1:3),Y(1:3)REAL KE(1:3,1:3)DO I=1,ND DO J=1,IW SK(I,J)=0 ENDDO ENDDO DO L=1,NE DO J=1,3 I=I3(J,L)X(J)=XY(1,I)Y(J)=XY(2,I)ENDDO CALL UKE1(X,Y,KE)DO 1 J=1,3 NJ=I3(J,L)DO 1 K=1,J NK=I3(K,L)IF(NJ.LT.NK)GOTO 11 NK=(NK-NJ+IW)SK(NJ,NK)=SK(NJ,NK)+KE(J,K)GOTO 1 11 NJ=(NJ-NK+IW)SK(NK,NJ)=SK(NK,NJ)+KE(J,K)NJ=NJ+NK-IW 1 CONTINUE ENDDO RETURN END !UKE1 计算单元系数矩阵 KE SUBROUTINE UKE1(X,Y,KE)REAL X(1:3),Y(1:3),C(1:3),B(1:3)REAL KE(1:3,1:3)C(1)=Y(2)-Y(3)C(2)=Y(3)-Y(1)C(3)=Y(1)-Y(2)B(1)=X(3)-X(2)B(2)=X(1)-X(3)B(3)=X(2)-X(1)S=2.*(C(1)*B(2)-C(2)*B(1)DO I=1,3 DO J=1,I KE(I,J)=(C(I)*C(J)+B(I)*B(J)/S ENDDO ENDDO RETURN END !UB1 加上第一类边界条件 SUBROUTINE UB1(ND1,NB1,U1,ND,IW,SK,U)REAL NB1(1:ND1),U1(1:ND1),SK(1:ND,1:IW),U(1:ND)DO I=1,ND U(I)=0.ENDDO DO I=1,ND1 J=NB1(I)SK(J,IW)=SK(J,IW)*1.E10 U(J)=SK(J,IW)*U1(I)ENDDO RETURN END !LDLT 解方程组 SUBROUTINE LDLT(A,N,IW,P,IE)REAL A(1:N,1:IW),P(1:N)DO 15 I=1,N IF(I.LE.IW)GO TO 20 IT=I-IW+1 GO TO 30 20 IT=1 30 K=I-1 IF(I.EQ.1)GO TO 40 DO 25 L=IT,K IL=L+IW-I B=A(I,IL)A(I,IL)=B/A(L,IW)P(I)=P(I)-A(I,IL)*P(L)MI=L+1 DO 25 J=MI,I IJ=J+IW-I JL=L+IW-J 25 A(I,IJ)=A(I,IJ)-A(J,JL)*B 40 IF(A(I,IW).EQ.0.)GO TO 100 15 CONTINUE DO 45 J=1,N IF(J.LE.IW)GOTO 60 IT=N-J+IW GOTO 70 60 IT=N 70 I=N-J+1 P(I)=P(I)/A(I,IW)IF(J.EQ.1)GOTO 45 K=I+1 DO 65 MJ=K,IT IJ=I-MJ+IW 65 P(I)=P(I)-P(MJ)*A(MJ,IJ)45 CONTINUE IE=0 GOTO 110 100 IE=1 110 RETURN END