有限元與邊界元(上)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)應(yīng)用地球物理系列課程有限元與邊界元教 案 李 貅 20068 教學(xué)基本情況課程學(xué)時課程總學(xué)時50，講授40學(xué)時，實(shí)驗(yàn)10課程要求 （1）掌握有限元法和邊界元法的基本理論、基本算法； （2）具有將實(shí)際地球物理問題中的一般邊值問題，利用有限元和邊界元的基本原理，轉(zhuǎn)化為有限元方程和邊界積分方程的能力； （3）了解有限元法和邊界元法的基本編程思想； （4）初步具有運(yùn)用有限元法和邊界元法解釋實(shí)際地球物理問題的能力。 課程的重點(diǎn)和難點(diǎn)（1）本課程的重點(diǎn)：有限元和邊界元的基本理論和基

2、本算法；解二維拉普拉斯方程的有限元法和邊界元法；解二維赫姆霍茲方程的有限元法和邊界元法；有限元法和邊界元法的應(yīng)用。（2）本課程的難點(diǎn)：各種邊值問題的變分原理；里茲伽略金法；基本解與格林公式。課程學(xué)時分配課程內(nèi)容學(xué) 時 數(shù)備 注總學(xué)時講授實(shí)驗(yàn)上機(jī)第一部分 有限元法28226第一章 有限元法數(shù)學(xué)基礎(chǔ)88第二章 有限元方法44第三章 解二維拉普拉斯方程的有限元法642第四章 解二維赫姆霍茲方程的有限元法1064第二部分 邊界元法22184第一章 邊界元法數(shù)學(xué)基礎(chǔ)66第二章 邊界單元法44第三章 解二維拉普拉斯方程的邊界元法642第四章 解二維赫姆霍茲方程的邊界元法642合 計504010第一部分 有

3、限元法序要求學(xué)生明確4個問題什么是有限元法_有限元法-是以變分原理和剖分插值為基礎(chǔ)的數(shù)值計算方法。有限元法實(shí)現(xiàn)過程的哲學(xué)思想邊值問題 變分問題（泛函極值問題） 剖分插值 （變分原理） （有限元法） 節(jié)點(diǎn)上未知量的高階線性方程組 求出節(jié)點(diǎn)場值。 有限元法的優(yōu)缺點(diǎn)優(yōu)點(diǎn)：1）于物性分布復(fù)雜的地球物理問題 2）解題過程規(guī)范缺點(diǎn)：全區(qū)域剖分單元和節(jié)點(diǎn)數(shù)目多現(xiàn)性代數(shù)方程組階數(shù)大對于無界區(qū)域問題須用大型計算機(jī)4、有限元法的應(yīng)用1）應(yīng)用條件-給出正確的邊值問題；有計算程序和計算機(jī)。2）應(yīng)用領(lǐng)域-地球物理中：位場延拓；重、磁、電、震、熱等正演。 工程應(yīng)用：3）解決地質(zhì)問題的特點(diǎn)-定量解釋4) 有限元法在地球物理

4、中所起的作用-解決了從前無法計算的地球物理問題；為地球物理反演奠定了基礎(chǔ)；提高了地球物理勘探的地質(zhì)效果。1.圍巖穩(wěn)定性分析圍巖分析中采用理想性模型模擬地下洞室開挖過程中洞壁的應(yīng)力和變形，并對錨固、襯石切等工程加固措施設(shè)計了專門的單元。該軟件采用Windows 風(fēng)格的用戶使用界面，工程技術(shù)人員只需輸入必要的工程結(jié)構(gòu)參數(shù)，就可很快地得到有限元計算結(jié)果。這軟件包對于洞室地址選擇，工程方案優(yōu)化有直接的指導(dǎo)意義。利用該軟件已應(yīng)用于北京十三陵水庫發(fā)電廠房的設(shè)計和云南大朝山水電站地下廠房的選址。2.小灣拱壩動力分析在有限元計算過程中，采用非連續(xù)體模型處理裂縫和斷層，考慮了壩體、基礎(chǔ)、壩肩的相互作用，考慮了水

5、載對計算的影響。該程序成功地模擬了小灣拱壩在地震作用的穩(wěn)定性。3.盆地模擬是利用FEPG系統(tǒng)開發(fā)的盆地演化有限元程序，已由中國海洋石油總公司應(yīng)用在我國南海的鶯瓊盆地和珠江口盆地、東海盆地、渤海灣盆地。應(yīng)用單位認(rèn)為該程序已成為石油地質(zhì)研究和目標(biāo)勘探的得力工具，為降低油氣勘探風(fēng)險發(fā)揮了重要作用。4.直流電極測井分析5.鉆孔測井分析6.塔里木盆地油氣運(yùn)移研究7地下廠房洞室群的三維圍巖穩(wěn)定分析 工程背景進(jìn)行索風(fēng)營電站的建設(shè)和進(jìn)一步優(yōu)化索風(fēng)營電站地下工程的設(shè)計和施工。通過模擬計算揭示出索風(fēng)營圍巖穩(wěn)定狀態(tài)，進(jìn)行圍巖穩(wěn)定評價，為設(shè)計的優(yōu)化和施工程序的安排提供指導(dǎo)。 地下廠房結(jié)構(gòu)圖 利用組合網(wǎng)格技術(shù)處理軟弱面

6、軟弱結(jié)構(gòu)面三維造型 8.三維電阻率測深有限元正演模擬中的邊界影響 二、三維低阻體模型 模型設(shè)計為區(qū)域的中心有一埋深1米，4m4m3m的低阻體，電阻率為50，為背景電阻率的1/2。邊界單元電阻率為2 和5000，分別為非邊界區(qū)單元電阻率的1/50和50倍，對比三種不同方向邊界單元電阻率變化的測深結(jié)果。 第一章 有限元法數(shù)學(xué)基礎(chǔ)變分法1.1泛函與變分問題1、泛函的概念：就是函數(shù)的函數(shù)。如：有一函數(shù) y=y(x),如果v又是y=y(x)的函數(shù)，則v=vy=vy(x)稱v為y的泛函。泛函和復(fù)合函數(shù)的區(qū)別：復(fù)合函數(shù)y=sin(x), z=y2 ; 泛函v=vy(x)2、泛函極值的概念變分問題以例子來說明

7、：例1：連接兩點(diǎn)弧長的最短線可分為兩個問題a)弧長問題泛函 A、B 兩點(diǎn)弧長 y=y（x）yxAB b）弧長最短線問題泛函極值問題求滿足下列條件：的y，這成為泛函的極值問題。例2：質(zhì)點(diǎn)沿曲線自由下滑的時間xy AVB質(zhì)點(diǎn)從A沿y（x）滑至B點(diǎn)，所需時間為 程t為y的泛函.A、B兩點(diǎn)的最速下降問題，即滿足下列條件：的y，這就是變分問題。1.2泛函極值與變分變分問題就是泛函的極值問題。泛函極值的計算方法類似于函數(shù)的極值的計算方法。變分的概念 在泛函中，自變量y（x）的增量是指滿足同一邊界的兩個之差成為自變量y的變分。 應(yīng)注意的兩點(diǎn)：變分與微分的區(qū)別 變分是對應(yīng)于同一個x的兩個之差 微分是x變化引起

8、的y的微分若固定，則有無限多種泛函的變分： 泛函的增量 泛函的變分泛函的極值： 變分的計算公式 泛函的極值條件 與函數(shù)的極值一樣，判別泛函的極大與極小還要看二階變分而定。1.3尤拉方程（變分原理）在以上各節(jié)中，我們有了變分問題的概念，但如何求變分問題是有限元法的關(guān)鍵。尤拉方程法是求解變分問題的方法之一，具體思路是： 首先將變分問題轉(zhuǎn)變?yōu)槲⒎址匠蹋ㄓ壤匠蹋?，然后解尤拉方程，得到變分問題的解。尤拉方程的建立將變分問題寫成一般形式： 設(shè)v在上取極值，任取一條與接近的曲線，y的變分為 考慮到：兩端點(diǎn)處的變分為零 變分對x的導(dǎo)數(shù)就是導(dǎo)數(shù)的變分 X1X2下面求泛函的變分：根據(jù)變分的計算公式有令 所以有

9、因?yàn)?所以用分部積分法計算上式中的第二項，由于，有在端點(diǎn)的變分為零，代入上式，得 所以 由泛函的極值條件，得 由的任意性，除兩點(diǎn)外，在（）內(nèi)，所以必有 -這就是尤拉方程。尤拉方程是泛函取極值時是必須滿足的方程，可通過積分的方法求得方程的通解y=y(x,c1,c2),由邊界條件可確定常數(shù)c1 , c2 。 尤拉方程展開式 由此可以看出，如果存在一個變分問題，可通過尤拉方程求出變分問題的解，到此為止我們解決了變分問題。但是如何將邊值問題轉(zhuǎn)化為變分問題還需要討論。 下面給出求變分問題的三個例子（見教材P10-13） 例一 求下列變分問題 例二 兩點(diǎn)的最短線問題 例三 最速下降問題極小位能原理和虛功原

10、理（兩個原理）從另一個角度出發(fā)建立起將邊值問題轉(zhuǎn)化為變分問題的方程。（1）弦的平衡問題兩點(diǎn)邊值問題 其中T為弦的張力（正常數(shù)），f(x)為外荷載垂直向下作用與弦。該問題的解就是顯得平衡位置。 由力學(xué)的“極小位能原理”可知，弦的平衡位置y*(x)是在滿足上述邊界條件的一切可能的位置y(x)中，使弦的總位能取最小者。弦處于某一位置y(x)時的總位能包括兩部分：一部分為弦的內(nèi)能，另一部分為外力所做的功。 總位能 其中稱為二次泛函。將上述兩個能量積分進(jìn)行變形，有 二次泛函又可寫成 其中（Ly,y）稱為Ly與y的內(nèi)積，,L稱為算子。上式對一切一維橢圓形方程的兩點(diǎn)邊值問題都適用。由于于是二次泛函可寫為 a

11、(y,y)稱為雙線形泛函，它具有雙線性； 2）對稱性； 3）正定性。為了確定弦的平衡位置，便導(dǎo)致兩個不同形式的數(shù)學(xué)問題：兩點(diǎn)邊值問題、變分問題，顯然二者之間應(yīng)具有某種等價關(guān)系，這就是一般要建立變分原理的最簡模型。（2）極小位能原理與虛功原理 定理1：若y*(x)是邊值問題 的解，則y*(x)使達(dá)到極小值；反之，若y*(x)使達(dá)到極小，則y*(x)是以上邊值問題的解。極小位能原理。定理2：若y*(x)是邊值問題 的解，則它必滿足虛功方程 反之，若y*(x)滿足虛功方程，則y*(x)是以上邊值問題的解。虛功原理。 虛功原理與極小位能原理相比，它不需要雙線性泛函a(y,z)具有對稱正定性，因此，虛功

12、原理的應(yīng)用范圍更廣。有了這兩個原理，我們可以將邊值問題轉(zhuǎn)換為變分問題，實(shí)際上虛功方程是尤拉方程的積分形式，兩者是一致的。1.4 依賴多個自變量的函數(shù)的泛函的變分問題設(shè)函數(shù)u是兩個自變量x,y的函數(shù)：u=u (x , y),且u（x , y）在區(qū)域的邊界上的值是給定的。求滿足這一邊界條件時，泛函 取極值的u (x , y)。 經(jīng)推導(dǎo)可知二維問題的尤拉方程為 解這個方程，得到帶積分常數(shù)的二維函數(shù)u（x ，y），帶入邊界條件后，解出積分常數(shù)，最后得到二維函數(shù)u（x ，y）。 例1：求變分問題相對應(yīng)的邊值問題（教材第15 頁） 設(shè)函數(shù)u是三個自變量x ，y，z的函數(shù)：u=u（x，y，z）在區(qū)域的邊界上

13、的值是給定的。求滿足這一邊界條件時，泛函經(jīng)推導(dǎo)可知三維問題的尤拉方程為 解這個方程，得到帶積分常數(shù)的三維函數(shù)u（x，y）, 帶入邊界條件后，解出積分常數(shù)，最后得到二維函數(shù) u（x，y）。 例2：求變積分問題相應(yīng)的邊值問題1.5 用里茲伽略金法解變分問題現(xiàn)介紹兩種古典的解變分問題方法。里茲法里茲法用一個線性獨(dú)立、完備的函數(shù)系 （i=1，-,n）,中的若干個函數(shù)的線性組合 作為泛函的試探解，其中，是待定系數(shù)，試探解應(yīng)在整個區(qū)域上有定義，而且滿足邊界條件。將試探解代入泛函中：，泛函v成為這些系數(shù)，i=1, -,n的函數(shù)：選擇使v取極值，v必須滿足 由上式得到含n個未知數(shù)的n個線性代數(shù)方程組，解方程得

14、，代入線性組合式，得近似解。伽略金法它是解虛功方程型的變分問題的一種近似方法，它的目標(biāo)也是求一組系數(shù)使 滿足虛功方程 由此虛功方程變?yōu)?考慮到j(luò)分別取1,2,n的每一項時,上式可分寫為 也可寫成 它和里茲法導(dǎo)出的代數(shù)方程組完全相同，因此習(xí)慣上成為里茲伽略金方程。以上方程為 解以上方程組，可確定。第二章 有限元方法方法選取基函數(shù)是比較困難的，而有限元法就在這一點(diǎn)上改進(jìn)和發(fā)展了方法，它用剖分和插值的辦法構(gòu)造基函數(shù)，下面以泊松方程的二維邊值問題為例說明有限元法的求解過程。2.1二維自然坐標(biāo)1、泊松方程的變分問題 邊值問題： 等價的變分問題：2、域剖分 用n個節(jié)點(diǎn)將區(qū)域剖分成m個單元，其剖分規(guī)則為：各

15、單元的節(jié)點(diǎn)只能與相鄰單元的節(jié)點(diǎn)相重，而不能成為相鄰單元的內(nèi)點(diǎn)。每個單元內(nèi)介質(zhì)的參數(shù)為常數(shù)。三角形單元的三條邊長之間差異不宜過大。在場變化劇烈的地方，應(yīng)使單元剖分的細(xì)一些。二維自然坐標(biāo)定義 設(shè)第e個單元所在平面的三個節(jié)點(diǎn)按逆時針方向編號，分別為：i，j，m，坐標(biāo)為對應(yīng)的節(jié)點(diǎn)函數(shù)值這三點(diǎn)組成三角形單元，其面積用表示，三角形中的點(diǎn)p（x，z）與i，j，m，三點(diǎn)的連線，將分割成三個小三角形其面積分別為點(diǎn)p在單元的位置，可表示為 稱為二維自然坐標(biāo)，其特點(diǎn)為定義i點(diǎn)：j點(diǎn)：m點(diǎn)：jipmxz均為線性函數(shù)。 其中 2)插值函數(shù) 假設(shè)單元內(nèi)u（x，z）是線性變化的u（x，z）=a+bx+cz a，b，c為系數(shù)

16、，三個節(jié)點(diǎn)上有由此得單元e內(nèi)函數(shù)u的線性插值表示式其中：稱為形函數(shù)。經(jīng)對比知形函數(shù)等于面坐標(biāo)故2.2 單元分析整個求解區(qū)域的積分vu可分解為單元積分veu 之和 取單元積分考慮到將泛函表達(dá)式寫成矩陣形式其中ni為頂點(diǎn)的單元個數(shù)?？傮w合成 把單元泛函累加求極值 根據(jù)極值理論 得 寫成矩陣形式 解方程就可以求得節(jié)點(diǎn)上的場值。2.3高斯數(shù)值積分在單元積分中，經(jīng)常遇到下列積分類型計算這類積分最常用的方法是高斯積方法，高斯積分法又稱最高代數(shù)精度求積公式。高斯數(shù)值積分的概念： 函數(shù)在區(qū)間-1，1的數(shù)值積分的一般形式 其中是區(qū)間-1，1中的幾個積分點(diǎn)，Ak是加權(quán)系數(shù)。如果給定時，選擇適當(dāng)?shù)腁k,，可是積分使

17、對任意的次多項式是精確的。 在給定后，上式中有n個待定系數(shù)Ak，令上式分別對1，是精確的，從而得到關(guān)于Ak的n個線性方程組，從方程組中可解出Ak來。數(shù)學(xué)上已證明，勒讓德（Legendre）多項式 的根作為積分點(diǎn)，而且求積分系數(shù)按下式計算 由此做出的高斯積分公式 高斯求積公式表見書P85。二維高斯積分公式三維高斯積分公式第三章 解二維拉普拉斯方程的有限元法31位場向上延拓與變換的有限元方法首先應(yīng)明確三個問題：什么是位場延拓根據(jù)地表實(shí)測的場值計算上部空間的場值，稱為向上延拓。在無源空間中，引力場、磁場都是位場，滿足拉普拉斯方程，所以，以上問題的延拓統(tǒng)稱為位場延拓。位場延拓的作用如果地下空間同時存在

18、大場源和小場源，則上部空間的場值較地表的場值更能壓制小場源的影響，因而，向上延拓有利于壓制干擾突出背景異常，向下延拓有利于突出局部異常。延拓的原理通過解邊值問題可以實(shí)現(xiàn)延拓；用格林公式可以實(shí)現(xiàn)延拓；利用惠根斯原理可以實(shí)現(xiàn)延拓。對于復(fù)雜界面的延拓需用有限元或邊界元法來實(shí)現(xiàn)。1、邊值問題1）方程與能量一般位函數(shù)u都滿足拉普拉斯方程引力場問題：引力場的垂直分量g時引力位u在垂直方向-y的負(fù)梯度磁場問題：磁場的垂直分量Z時磁位u在垂直方向-y的負(fù)梯度磁場的水平分量H時磁位u在水平方向x的負(fù)梯度磁場T是磁位u在地磁場方向t的負(fù)梯度2）邊界條件Sn 上半空間的邊界有一個無窮大邊界，與地表邊界組成封閉的邊界

19、。但上的u是未知的，可用以下幾種方法處理： SHAPE * MERGEFORMAT 將取得足夠遠(yuǎn)，近似地認(rèn)為，這樣會使節(jié)點(diǎn)數(shù)增加計算工作量增大。將近似看作水平線，用向上延拓公式，由上的u計算出上的u。這樣可以減輕計算量，但是近似地。以上兩種方法，得到第一類邊界條件其中是已知函數(shù)。近似認(rèn)為地表的u是地下某點(diǎn)處的集中場源s產(chǎn)生的，上u的分布規(guī)律為其偏導(dǎo)數(shù)利用以上兩式消去常數(shù)后得這就是第三類邊界條件。另一種情況，如果在u的分布規(guī)律為則第三類邊界條件為以上的邊值問題歸結(jié)為2、變分問題第一類邊界條件的邊值問題對應(yīng)的變分問題為第三類邊值條件的邊值問題其中為第二類邊界條件，為第三類邊界條件。對應(yīng)的變分問題由

20、于在泛函式取極值過程中，邊界條件已考慮在內(nèi)，不需解方程組時再考慮，這種第二類、第三類邊界條件成為自然邊界條件?；旌线吔鐥l件的邊值問題對應(yīng)的變分問題3、有限單元法（1）區(qū)域劃分用三角形元對整個區(qū)域進(jìn)行剖分，剖分后對節(jié)點(diǎn)和單元進(jìn)行編號，將節(jié)點(diǎn)的x,y坐標(biāo)和單元的節(jié)點(diǎn)號列表xy12340123456789101112131415161718192021（1）（2）（3）（4）（5）（6）（7）（8）（9）（10）（11）（12）（13）（14）（15）（16）（17）（18）（20）（21）（22）（23）（24）（25）（26）（27）（28）（19）節(jié)點(diǎn)的坐標(biāo)節(jié)點(diǎn)號12NDxX1X2XNDyY1

21、Y2YND單元的節(jié)點(diǎn)號單元號12NEii1i2iNEjj1j2jNEkk1k2kNE節(jié)點(diǎn)的場值序號12ND1場值u1u2uNE1(2)線性插值假設(shè)單元內(nèi)u(x,z)是線性變化的u(x,z)=a+bx+cz a,b,c為系數(shù)，其中：（3）單元分析整個求解區(qū)域的積分Fu可分解為單元積分Feu之和考慮到將泛函表達(dá)式寫成矩陣形式 ， 其中ni為以I為頂點(diǎn)的單元個數(shù)。(4)總體合成把單元泛函累加其中是總體系數(shù)矩陣。（5）求變分對上式求變分 其中，除第一類邊界條件的節(jié)點(diǎn)的外，其余節(jié)點(diǎn)上的，因?yàn)镵是對稱矩陣，有所以由于，所以解方程就可求得節(jié)點(diǎn)上的場值。（6）解方程有限元法程序框圖輸入原始數(shù)據(jù)節(jié)點(diǎn)總數(shù)ND，單

22、元總數(shù)NE；節(jié)點(diǎn)坐標(biāo)XY(2,ND)；單元編號I3(3,NE)邊界節(jié)點(diǎn)號NB1(ND1)；邊界節(jié)點(diǎn)電場值U1(ND1)計算半帶寬（調(diào)用子程序MBW）代入第一類邊界條件（調(diào)用子函數(shù)UB1）形成總體系數(shù)矩陣（調(diào)用子程序UK1）解線性方程組（調(diào)用子程序LDLT）輸出數(shù)字解PROGRAM maincharacter *20filename1PARAMETER(ND=15,NE=16,ND1=12) DIMENSION I3(3,16),XY(2,15),U1(12),NB1(12),SK(15,15),U(15)OPEN(11,FILE=NE.txt,STATUS=old)DO 1 J=1,NE1RE

23、AD(11,*)(I3(I,J),I=1,3)READ(11,*)(U1(i),i=1,ND1)DO 2 J=1,ND2READ(11,*)(XY(I,J),I=1,2)READ(11,*)(NB1(I),I=1,ND1)CLOSE(11)CALL MBW(NE,I3,IW)CALL UK1(ND,NE,IW,I3,XY,SK)CALL UB1(ND1,NB1,U1,ND,IW,SK,U) CALL LDLT(SK,ND,IW,U,IE)write(*,*)請輸入保存計算后的數(shù)據(jù)文件名：read(*,*)filename1open(6,file=filename1,status=unknown

24、) write(6,510)Uclose(6)510format(/3e15.6)ENDSUBROUTINE MBW(NE,I3,IW)DIMENSION I3(3,NE)IW=0DO 10 I=1,NEM=MAX(IABS(I3(1,I)-I3(2,I),IABS(I3(2,I)-I3(3,I), &IABS(I3(3,I)-I3(1,I)IF(M+1.GT.IW)IW=M+110CONTINUERETURNENDSUBROUTINE UK1(ND,NE,IW,I3,XY,SK)DIMENSION I3(3,NE),XY(2,ND),SK(ND,IW)DIMENSION X(3),Y(3)R

25、EAL KE(3,3)DO 10 I=1,NDDO 10 J=1,IW10SK(I,J)=0.DO 20 L=1,NEDO 30 J=1,3I=I3(J,L)X(J)=XY(1,I)30Y(J)=XY(2,I)CALL UKE1(X,Y,KE)DO 40 J=1,3NJ=I3(J,L)DO 40 K=1,JNK=I3(K,L)IF(NJ.LT.NK)GOTO 50NK=NK-NJ+IWSK(NJ,NK)=SK(NJ,NK)+KE(J,K)GOTO 4050NJ=NJ-NK+IWSK(NK,NJ)=SK(NK,NJ)+KE(J,K)NJ=NJ+NK-IW40CONTINUE20CONTINUER

26、ETURNENDSUBROUTINE UKE1(X,Y,KE)DIMENSION X(3),Y(3),A(3),B(3)REAL KE(3,3)A(1)=Y(2)-Y(3)A(2)=Y(3)-Y(1)A(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.*(A(1)*B(2)-A(2)*B(1)DO 10 I=1,3DO 10 J=1,I10KE(I,J)=(A(I)*A(J)+B(I)*B(J)/SRETURNENDSUBROUTINE UB1(ND1,NB1,U1,ND,IW,SK,U)DIMENSION NB1(ND1),

27、U1(ND1),SK(ND,IW),U(ND)DO 10 I=1,ND10U(I)=0.DO 20 I=1,ND1J=NB1(I)SK(J,IW)=SK(J,IW)*1.E1020U(J)=SK(J,IW)*U1(I)RETURNENDSUBROUTINE LDLT(A,N,IW,P,IE)DIMENSION A(N,IW),P(N)DO 15 I=1,NIF(I.LE.IW) GOTO 20IT=I-IW+1GOTO 3020IT=130K=I-1IF(I.EQ.1)GOTO 40DO 25 L=IT,KIL=L+IW-IB=A(I,IL)A(I,IL)=B/A(L,IW)P(I)=P(I)

28、-A(I,IL)*P(L)MI=L+1DO 25 J=MI,IIJ=J+IW-IJL=L+IW-J25A(I,IJ)=A(I,IJ)-A(J,JL)*B40IF(A(I,IW).EQ.0.)GOTO 10015CONTINUEDO 45 J=1,NIF(J.LE.IW)GOTO 60IT=N-J+IWGOTO 7060IT=N70I=N-J+1P(I)=P(I)/A(I,IW)IF(J.EQ.1)GOTO 45K=I+1DO 65 MJ=K,ITIJ=I-MJ+IW65P(I)=P(I)-P(MJ)*A(MJ,IJ)45CONTINUEIE=0GOTO 110100IE=1110RETURNE

29、ND32二維均勻電場的電阻率法的有限元算法1、邊值問題假定地下構(gòu)造是二維的，垂直構(gòu)造走向做一剖面，選取足夠大的矩形區(qū)域?yàn)檠芯繀^(qū)域，AD是水平地面，AB，BC，CD代表無窮遠(yuǎn)邊界區(qū)域中的初始電場是水平均勻場，其強(qiáng)度為E0，異常電場在無窮遠(yuǎn)邊界上為零，電位的邊值問題為 CAAAAAAAADB2、變分問題根據(jù)極小位能原理以上方程對應(yīng)的泛函由多元函數(shù)的變分求和法知：若泛函對于上式故利用所以利用場論中的關(guān)系式可推得取故以上變分表達(dá)式可寫為 其中第二、四項積分為零（由方程），第一、三項積分可變成邊界積分。根據(jù)交接面條件知于是可見，由于內(nèi)邊界條件是自然條件，因此在變分中不出現(xiàn)，變分只與外部邊界有關(guān)。將上式可

30、寫成移項可得 因此，電位的邊值問題與下列變分問題等價3、有限元法（1）矩形單元、雙線性差值區(qū)域劃分：用規(guī)則的矩形網(wǎng)格對區(qū)域進(jìn)行剖分，節(jié)點(diǎn)編號和單元編號如圖9516210073118412雙線性插值：y2314（0，1）2（1，0）2(a)2314(b)xab是母單元，(b)是子單元。兩個單元間的坐標(biāo)變換關(guān)系其中是子單元中點(diǎn)坐標(biāo)，a,b是子單元的兩個邊長。微分關(guān)系為構(gòu)造如下形函數(shù)其中，是點(diǎn)i(i=1,2,3,4)的坐標(biāo)，形函數(shù)滿足其中j代表點(diǎn)號。單元中u的插值函數(shù)是其中(i=1,2,3,4)是單元四個頂點(diǎn)的待定值。這顯然是關(guān)于x,y的線性函數(shù)。單元分析：單元上的積分其中考慮到所以積分得其中 經(jīng)計

31、算可得現(xiàn)在考慮邊界積分部分，若單元的一個邊落在AB上，則邊界積分其中同理CD邊界上的積分其中單元積分可寫出總體集成：將全部單元相加得求變分：由于，所以這是含有ND個單元的ND個方程的線性代數(shù)方程組。第四章 解二維赫姆霍茲方程的有限元法赫姆霍茲方程是地球物理中常見的方程，比如穩(wěn)定電流場中的點(diǎn)源二維問題，在這種情況下，為了求解方便，往往延拓對稱軸方向進(jìn)行付氏變換，進(jìn)而將拉普拉斯方程轉(zhuǎn)化為赫姆霍茲方程。另外，頻率域的電磁場問題均滿足赫姆霍茲方程。4.1點(diǎn)源二維電場的計算方法點(diǎn)源是指點(diǎn)電源供電，產(chǎn)生的電場是三維的。二維是有一定走向的二維構(gòu)造其電性特征是二維分布?？傮w來講，點(diǎn)源二維構(gòu)造的電場實(shí)質(zhì)上是三維

32、的。對于這樣的問題可以用付氏變換的方法，將三維電場變換成帶參數(shù)的二維問題來處理，使得計算更方便。1、邊值問題設(shè)在地面A點(diǎn)，置一電流強(qiáng)度為I的點(diǎn)電源，地下構(gòu)造成二維分布，如圖所示。yxzrA對u在z方向上進(jìn)行付氏變換：由付氏變換的微分性質(zhì) 變換后的二維含參邊值問題 （第三類邊界條件）2、變分問題首先構(gòu)造一個泛函（注：解釋應(yīng)用極小位能原理和虛功原理來構(gòu)造泛函）下一步要用變分原理驗(yàn)證我們構(gòu)造的變分問題與邊值問題是等價的，這樣就可以建立變分方程。求變分：由方程知將上式中的外部邊界條件代入所以所以二維邊值問題與下列變分問題等價3、有限單元法用有限單元法求解點(diǎn)源二維電場的變分方程所用的區(qū)域剖分、插值等方法與求解二維均勻電場變分方程所用的方法相同。（略）4、傅立葉反變換用以上方法解得的各節(jié)點(diǎn)的U，該U對應(yīng)于特定的波數(shù)k。如果計算了對應(yīng)于不同k的一組U可用傅立葉反變換計算三維空間的電位u。關(guān)鍵問題：數(shù)值計算方法