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

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

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

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

4、的單元和節(jié) 點(diǎn)數(shù)目多，最后得到的線性代數(shù)方程組很大。特別是三維問(wèn)題和地球物理中常遇 到的無(wú)解區(qū)域問(wèn)題，需要中、大型計(jì)算機(jī)，才能完成有限單元法的計(jì)算。這是有 限單元法的主要缺點(diǎn)。 基本步驟 用三角單元對(duì)區(qū)域進(jìn)行剖分： 用三角單元對(duì)整個(gè)區(qū)域進(jìn)行剖分，因?yàn)槿菃卧倪吶菀讛M合地形線的形 狀。三角形的頂點(diǎn)稱(chēng)為節(jié)點(diǎn)，用節(jié)點(diǎn)上的離散場(chǎng)值來(lái)近似場(chǎng)值的連續(xù)分布。剖分 后對(duì)單元和節(jié)點(diǎn)進(jìn)行編號(hào)，次序號(hào)(i，j， m)按逆時(shí)針?lè)较蚺帕?。第一?lèi)邊界 條件的節(jié)點(diǎn)上的場(chǎng)值是已知的，其余節(jié)點(diǎn)上的場(chǎng)值是待求的； 線性插值： 在三角單元內(nèi)假定場(chǎng)值 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是對(duì)稱(chēng)矩陣。 4 八 總體合成： 將單元的場(chǎng)值列向量 Ue擴(kuò)展成全體節(jié)點(diǎn)的場(chǎng)值向量 U， U=(u1u2Uujumund)，按照節(jié)點(diǎn)的總體序號(hào)，將單元系數(shù)矩陣中的各元放在 ke的相應(yīng)行與列的交叉位置上，其余位置的元為零，這樣單元積分可寫(xiě)成 1 T 1 T Fe(U)：二 Ueke% U 2 2 F(u)= )2dxdy , 各單元積分相加時(shí)，只要將 ke相加即可； 求變分： 通過(guò)以上四個(gè)步驟，已經(jīng)將連續(xù)函數(shù) g 的泛函，離散成各節(jié)點(diǎn) g 值的多 1 元函數(shù)：F( u)=

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

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

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

10、證 了有限單元法的可效性。 3、剖分后的節(jié)點(diǎn)分布及單元編號(hào)見(jiàn)下表： （1） 節(jié)點(diǎn)總數(shù) ND=33 ； （2） 單元總數(shù) NE=40; （3） 單元節(jié)點(diǎn)編號(hào) 表 1 單元 序號(hào) 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 測(cè)線 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é)點(diǎn)坐標(biāo) 表 2 節(jié)點(diǎn) 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) 第一類(lèi)邊界節(jié)點(diǎn)數(shù) ND1=24; (6) 第一類(lèi)邊界節(jié)點(diǎn)號(hào)和場(chǎng)值： 表 3 邊界節(jié) 點(diǎn)序號(hào) 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 點(diǎn)場(chǎng)值 (*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ù)書(shū)中給出的第一類(lèi)邊界條件、三角元剖分、線性插值的位場(chǎng)延拓得有限 單元程序框圖： 編制計(jì)算程序，并將已知的邊界節(jié)點(diǎn)數(shù)據(jù)輸入，來(lái)計(jì)算區(qū)域內(nèi)部的節(jié)點(diǎn)值（5、 8、11、14、17、20、23、26、29 點(diǎn)）。運(yùn)行程序后，得計(jì)算結(jié)果如下： 內(nèi)部節(jié)點(diǎn)號(hào) 數(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é) 通過(guò)計(jì)算出的解與用解析公式得到的解進(jìn)行比較， 誤差相差不大，證明了有 限元方法是一種有效的正演方法。附錄：計(jì)算機(jī)程序 地形及圖形剖分程序： 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; % 觀測(cè)點(diǎn) % %繪圖 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é)點(diǎn)坐標(biāo)值計(jì)算及比較的 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 !萬(wàn)有引力常量 X 仁 200.;Y1=-100. !球心 S 坐標(biāo)(x1,y2) R=10. !半徑 P1=1. !密度 DO 1=1,11 XX(1,I)=40.*(I-1) !觀測(cè)點(diǎn)坐標(biāo) xx YY(1,I)=20.*sin(0.02*XX(1

23、,l)+30. !觀測(cè)點(diǎn)坐標(biāo) yy XX(2,I)=XX(1,I) !網(wǎng)格點(diǎn)坐標(biāo) xx1 YY(2,l)=(300.+YY(1,l)/2. !網(wǎng)格點(diǎn)坐標(biāo) yy1 XX(3,I)=XX(1,I) !網(wǎng)格點(diǎn)坐標(biāo) xx2 YY(3,I)=300. !網(wǎng)格點(diǎn)坐標(biāo) yy2 ENDDO !計(jì)算異常場(chǎng)值 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é)點(diǎn)編號(hào) 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é)點(diǎn)號(hào)，”數(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 計(jì)算總體系數(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 計(jì)算單元系數(shù)矩陣 KE SUBROUTINE UKE1(X,Y ,KE)