地球物理中的有限單元法分解

1、地球物理算法技術(shù)（論文） 地球物理中的有限單元法 院 系:地球物理與信息技術(shù)院 姓 名:劉雅寧 學 號:2010120053 任課老師:張貴賓地球物理中的有限單元法 一、有限單元法的介紹 在地球物理理論計算中，存在著兩類基本問題：正問題和反問題。給定場源 的分布，求解場值的大小，這是正問題，或者稱為正演問題。 地球物理正演的數(shù)值計算方法，種類很多，最常用的有：有限差分法和有限 單元法。 有限單元法是 50 年代首先在彈性力學中發(fā)展起來的方法。主要優(yōu)點是，適 用于物性參數(shù)復雜分布的區(qū)域，但計算量大。隨著計算機技術(shù)的發(fā)展，有限單元 法在解決各個工程領(lǐng)域的許多數(shù)學物理問題中， 得到了廣泛的應(yīng)用，稱為

2、一種高 效、通用的計算方法。地球物理中的一些邊值問題，也采用了有限單元法，解決 了許多從前無法計算的地球物理問題。 有限單元法解決數(shù)學物理邊值問題的基本思路和過程如下： 1、 給出地球物理邊值問題中的偏微分方程和邊界條件（及初始條件） 。這一 點看起來似乎容易，但做起來并不容易，特別是邊界條件的給定。只有對地球物 理方法的原理和問題有深入的理解，才能給邊值問題中的偏微分方程和邊界條件 以正確的描述。 2、 將地球物理邊值問題轉(zhuǎn)變?yōu)橛邢拊匠?。實現(xiàn)這種轉(zhuǎn)變的主要數(shù)學工具是 變分法，用變分法得到的有限元法方程稱為泛函極值問題。 3、 用優(yōu)先單元法解決泛函極值問題其步驟大致如下： 把研究區(qū)域剖分成有

3、限個小單元，在每個單元上，把函數(shù)簡化成線性函數(shù)、 二次函數(shù)或高次函數(shù)，這稱為單元上函數(shù)的插值。用簡化后的函數(shù)計算每個單元 上的泛函。各單元之間，通過單元間節(jié)點上的函數(shù)值相互聯(lián)系起來。 對各單元的 泛函求和，獲得整個區(qū)域上的泛函。這樣，有限單元法將連續(xù)函數(shù)的泛函，離散 成各單元節(jié)點上函數(shù)值得泛函。根據(jù)泛函取極值的條件，得到各節(jié)點的函數(shù)值應(yīng) 滿足的線性代數(shù)方程組。解代數(shù)方程組，得到各節(jié)點的函數(shù)值。 有限單元法的主要優(yōu)點是，適用于物性復雜分布的地球物理問題， 而且，其 解題過程也比較規(guī)范化。這些優(yōu)點是有限單元法在地球物理中獲得廣泛的應(yīng)用。 但是，有限單元法是區(qū)域性方法，必須在全區(qū)域中進行剖分。剖分后

4、的單元和節(jié) 點數(shù)目多，最后得到的線性代數(shù)方程組很大。特別是三維問題和地球物理中常遇 到的無解區(qū)域問題，需要中、大型計算機，才能完成有限單元法的計算。這是有 限單元法的主要缺點。 基本步驟 用三角單元對區(qū)域進行剖分： 用三角單元對整個區(qū)域進行剖分，因為三角單元的邊容易擬合地形線的形 狀。三角形的頂點稱為節(jié)點，用節(jié)點上的離散場值來近似場值的連續(xù)分布。剖分 后對單元和節(jié)點進行編號，次序號(i，j， m)按逆時針方向排列。第一類邊界 條件的節(jié)點上的場值是已知的，其余節(jié)點上的場值是待求的； 線性插值： 在三角單元內(nèi)假定場值 u 是線性變化的，則 u=ax+by+c，u 還可表示為 u=Niui+叫Uj+

5、Nmum 其中 NjNjNm是形函數(shù)， 1 1 Ni= -(aiX+biy+Ci)，Nj= -(ajX+bjy+Cj)， 2L 2L 1 Nm= ( amX+bmy+Cm)， 2 也 ai =yj-ym，bi=Xm-Xj，Ci =Xjym-Xmyj， aj=ym-yi，bj=Xi-Xm，Cj=Xmyi -Xiym， am=yi-yj，bm=Xj-Xi，Cm=Xiyj-Xjyi， 1 -= -(aibj-ajbi)是三角元面積 2 單元分析： 將全區(qū)域的積分分解為單元 e 的積分之和 z)2dE hurdjJE (空)2+ 2 1 e 2 1 e 2 X 繼續(xù)推得單元 e 上的積分 其中u： =

6、(u u u m是單元的 u 值列向量；ke是單元系數(shù)矩陣， 1 kst= asat+bsbt)，s，t=i，j，m， ke是對稱矩陣。 4 八 總體合成： 將單元的場值列向量 Ue擴展成全體節(jié)點的場值向量 U， U=(u1u2Uujumund)，按照節(jié)點的總體序號，將單元系數(shù)矩陣中的各元放在 ke的相應(yīng)行與列的交叉位置上，其余位置的元為零，這樣單元積分可寫成 1 T 1 T Fe(U)：二 Ueke% U 2 2 F(u)= )2dxdy , 各單元積分相加時，只要將 ke相加即可； 求變分： 通過以上四個步驟，已經(jīng)將連續(xù)函數(shù) g 的泛函，離散成各節(jié)點 g 值的多 1 元函數(shù)：F( u)=

7、uTku，泛函的極值等于多元函數(shù)的極值，用多元函數(shù)求極值 2 的方法，對上式求微分， 1 1 F(u) : d (-uTku) = -(duTku+uTkdu)=0 2 2 因為 K 是對稱矩陣，有 duTku=uT kdu,所以:F (u du ku=0，由于 duT = 0，所以由上式得ku=0。得到含有 ND 個元的 ND 個方程聯(lián)立的線性 代數(shù)方程組； 解線性代數(shù)方程組： 首先將第一類邊界條件代入，通過定帶寬儲存的對稱帶型線性方程組， 解方程組，得到各節(jié)點的 u 值，至此有限單元法的求解過程結(jié)束。 三、實現(xiàn)過程 為了驗證有限單元法的有效性，我們設(shè)計一個規(guī)則形狀的地下礦體，給出模 型：

8、1 模型 密度均勻的水平半無限空間，一個均勻球體 S,球體半徑 R=10m，剩余密度 (T =1g/cm3,球心坐標(a，b) = (200, -100)。對于均勻球體來說，它與將其全 部剩余質(zhì)量集中在球心處的點的質(zhì)量所引起的異常完全一樣。若球心的埋藏深度 為 D，球的半徑為 R，剩余密度為，貝尼的剩余質(zhì)量為 M=(4 n R3C )/3,通過原 點的任意水平剖面上則重力異常的解析解表達式為： . GMD Q(x2 D2)3/2 設(shè)測線長 400m，高程變化(-200,300 )，地形設(shè)為一曲線： y=20*sin(0.02*x)+30，其中 x 為測線離原點的水平距離，y 為高程，則 S 引

9、起的 重力異常為： 4 二 R3G(y-b) g 3(x-a)2 (y-b)2)3/2 2、剖分 通過 Matlab 建模，得到地形曲線圖，再用三角單元對劃出的區(qū)域進行剖分 （程序附后），剖分后對單元和節(jié)點進行編號，并將節(jié)點的 xy 坐標和單元節(jié)點號 列表（表 1 和表 3），分別放在 XY（2，ND ）和 13（3，NE）兩個二維數(shù)組中。 剖分圖如下： -100 圖 1:剖分圖 根據(jù)重力異常的解析解表達式算出區(qū)域邊界上及區(qū)域內(nèi)部的節(jié)點值。 編制計算程 序，根據(jù)邊界上節(jié)點的場值，算出區(qū)域內(nèi)部節(jié)點的場值，將計算出來的內(nèi)部節(jié)點 場值與解析解比較，在有效數(shù)字內(nèi)，若計算值與理論值一致或相差不大， 則驗

10、證 了有限單元法的可效性。 3、剖分后的節(jié)點分布及單元編號見下表： （1） 節(jié)點總數(shù) ND=33 ； （2） 單元總數(shù) NE=40; （3） 單元節(jié)點編號 表 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 3 V J2 c15 . 18 廠1 少 0 亠3Q J -2OOo 50 100 150 200 測線 250 300 350

11、400 200 100 0 300 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

12、26 27 31 28 29 30 (4) 節(jié)點坐標 表 2 節(jié)點 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 172.1 300 49.9 174.9 300 43.5 171.7 4712 7355 9147 9573 0927 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 164.4 300 14.8 157.4 300 10.0 155.03 3

13、00 17.3 3252 1626 6395 3198 7671 836 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 300 32.3 166.1 300 45.8 172.9 300 49.7 174.8 300 8733 3098 6550 7336 3668 8717 9359 (5) 第一類邊界節(jié)點數(shù) ND1=24; (6) 第一類邊界節(jié)點號和場值： 表 3 邊界節(jié) 點序號 1 2 3 4 6 7 9 10 12 13 15 邊界節(jié) 2.676 2.023

14、 1.249 4.031 1.398 5.914 1.534 9.042 1.646 14.66 1.720 點場值 (*103) 8204 8112 8540 5208 0970 4786 9158 7687 9267 69839 8469 16 18 19 21 22 24 25 27 28 30 31 32 33 21.18 1.746 19.14 1.720 11.44 1.646 6.487 1.534 4.016 1.398 2.683 1.955 1.249 2479 7240 9463 8469 5633 9267 6173 9158 5414 0970 2691 5145 8

15、540 根據(jù)書中給出的第一類邊界條件、三角元剖分、線性插值的位場延拓得有限 單元程序框圖： 編制計算程序，并將已知的邊界節(jié)點數(shù)據(jù)輸入，來計算區(qū)域內(nèi)部的節(jié)點值（5、 8、11、14、17、20、23、26、29 點）。運行程序后，得計算結(jié)果如下： 內(nèi)部節(jié)點號 數(shù)值解 解析解 誤差 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.001

16、5380944 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 小結(jié) 通過計算出的解與用解析公式得到的解進行比較， 誤差相差不大，證明了有 限元方法是一種有效的正演方法。附錄：計算機程序 地形及圖形剖分程序： clc,clea

17、r x=0:400;y=-200:0; X Y=meshgrid(x,y); y 仁 20*sin(0.02*x)+30; % 地形 % Model=50*exp(-(200-X)A2/120-(-100-Y)A2/120); % 模型 % x 仁 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,li nestyle, non e) hold on contourf(X,Y ,Mode

18、l,400,linestyle,none) set(gca,tickle ngth,0.0 0.0,Fo ntName,Verda na,Fo ntWeight,Bold,fo ntsize,12)% colormap jet(32) for j=1:le ngth(xl) 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

19、i=1:3 hold on plot(x1,Y1(i,:),ko,x1,Y1(i,:),k,l in ewidth,1.2); end for j=1:le ngth(xl) hold on plot(zeros(3,1)+x1(j),Y1(:,j),ko,zeros(3,1)+x1(j),Y1(:,j),k,li newidth,1.2); end for i=1:le ngth(x1)-2 hold on plot(x1(i:i+1),Y1(1,i),Y1(2,i+1),k,li newidth,1.2); hold on plot(x1(i:i+1),Y1(3,i),Y1(2,i+1),

20、k,li newidth,1.2); end hold on plot(x1(le ngth(x1)-1:le ngth(x1),Y1(1,le ngth(x1)-1),Y1(2,le ngth(x1),k,li newidth,1.2); hold on plot(x1(le ngth(x1)-1:le ngth(x1),Y1(3,le ngth(x1)-1),Y1(2,le ngth(x1),k,li newidth,1.2); k=0; for j=1:le ngth(x1) for i=3:-1:1 k=k+1; text(X1(i,j)+2,Y1(i,j)+10, num2str(k

21、),color,r,fo ntsize,10,fo ntn ame, 華文中宋 ); end end axis(0 400 -200 400) 節(jié)點坐標值計算及比較的 Fortran 程序： PROGRAM seco nd 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(:,:),

22、A(:,:),I3(:,:) REAL,ALLOCA TABLE:P(:),NB1(:),U1(:),U(:) ND=33 NE=40 ND1=24 ALLOCA TE(U(1:ND),P(1:ND),NB1(1:ND1),U1(1:ND1),I3(1:3,1:NE) !地質(zhì)體的模型參數(shù) Pl=3.14159 !PI 常量 G=6.672E-11 !萬有引力常量 X 仁 200.;Y1=-100. !球心 S 坐標(x1,y2) R=10. !半徑 P1=1. !密度 DO 1=1,11 XX(1,I)=40.*(I-1) !觀測點坐標 xx YY(1,I)=20.*sin(0.02*XX(1

23、,l)+30. !觀測點坐標 yy XX(2,I)=XX(1,I) !網(wǎng)格點坐標 xx1 YY(2,l)=(300.+YY(1,l)/2. !網(wǎng)格點坐標 yy1 XX(3,I)=XX(1,I) !網(wǎng)格點坐標 xx2 YY(3,I)=300. !網(wǎng)格點坐標 yy2 ENDDO !計算異常場值 DO 1=1,3 DO J=1,11 ! 異常體 S 異常值 DG(I,J)=4.*PI*R*3.*P1*G*(YY(I,J)-Y1) & /(3.*(XX(l,J)-X1)*2.+(YY(l,J)-Y1)*2.)*(3.0/2.0)*1.0E9 WRITE(*,*) DG(I,j) ENDDO E

24、NDDO ! INPUT DA TA DO 1=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 ! 輸入?yún)?shù) I3，存放單元節(jié)點編號 OPEN(2,FILE=l3.txt) READ (2,*) DO I=1,NE

25、READ (2,*) I3(1,I),I3(2,I),I3(3,I) ENDDO CLOSE(2) WRITE (*,*) SUCCESS: In put 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 函數(shù) ! Call MBW CALL MBW(NE,I3,IW) WRITE (*,*) SUCCESS: Call MBW ALLOCA

26、TE(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

27、,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)節(jié)點號，”數(shù)值解，”解析解

28、，”誤差” DO I=1,ND WRITE (4,(l15,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 sec ond ! SUBROUTINE ! MBW 計算總體系數(shù)矩陣的半帶寬 SUBROUTINE MBW(NE,I3,IW) INTEGER M REAL I3(1:3,1:NE) IW=O DO I=1,NE M

29、=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 計算單元系數(shù)矩陣 KE SUBROUTINE UKE1(X,Y ,KE)