有限單元法的基本知識和地震波傳播正演模擬的應(yīng)用課件

1、地球物理數(shù)值計算方法地球物理數(shù)值計算方法Numerical Methods in Geophysics王彥賓 地球物理學(xué)系2008-2009學(xué)年第二學(xué)期1地球物理數(shù)值計算方法第六章 有限單元方法2有限單元法 地球內(nèi)部介質(zhì)，尤其是淺部，存在橫向非均勻結(jié)構(gòu)，包括分層、不規(guī)則形狀的塊體。由于解析方法不能給出這類復(fù)雜模型中的解，因此需要近似的數(shù)值求解方法。有限單元方法(Finite Element Method，F(xiàn)EM)是地球物理數(shù)值方法中另一種常用到的方法。 在有限元方法中，把計算域離散剖分為有限個互不重疊且相互連接的單元，在每個單元內(nèi)選擇基函數(shù)，用單元基函數(shù)的線性組合來逼近單元中的真解，整個計算域

2、上總體的基函數(shù)可以看為由每個單元基函數(shù)組成的，則整個計算域內(nèi)的解可以看作是由所有單元上的近似解構(gòu)成。根據(jù)所采用的權(quán)函數(shù)和插值函數(shù)的不同，有限元方法也分為多種計算格式。從計算單元網(wǎng)格的形狀來劃分，有三角形網(wǎng)格、四邊形網(wǎng)格和多邊形網(wǎng)格。 有限元方法最早應(yīng)用于結(jié)構(gòu)力學(xué)，后來隨著計算機(jī)的發(fā)展慢慢用于流體力學(xué)的數(shù)值模擬。 在地球物理領(lǐng)域，有限單元法應(yīng)用于地震波場模擬、地球動力學(xué)模擬等。由于網(wǎng)格劃分的靈活性，特別適用于非常復(fù)雜的模型，如自由地表、復(fù)雜邊界模型。3有限單元法簡介有限單元方法是將偏微分方程描述的連續(xù)問題進(jìn)行離散求解的一種數(shù)值方法。其基本原理是：用簡單的塊體構(gòu)造復(fù)雜的對象，或?qū)⒁粋€復(fù)雜的對象分為

3、用以處理的小塊體。例子：近似圓的面積 一個三角形的面積： 圓的面積： N為三角形的個數(shù)，當(dāng)N時，4有限單元法簡介有限單元方法包括以下基本步驟：1、根據(jù)變分原理或方程余量與權(quán)函數(shù)正交化原理，建立與微分 方程初邊值問題等價的積分表達(dá)式，這是有限元法的出發(fā) 點(diǎn)。2、將研究對象剖分為簡單的塊體，即單元。3、描述每一個單元中的物理量，確定單元基函數(shù)，將各個單元 中的求解函數(shù)用單元基函數(shù)的線性組合表達(dá)式進(jìn)行逼近。4、將所有單元通過節(jié)點(diǎn)組合到一起，按一定法則進(jìn)行累加，形 成總體有限元方程。5、處理邊界條件。6、采用適當(dāng)?shù)臄?shù)值計算方法求解根據(jù)邊界條件修正的總體有限 元方程組，可求得各節(jié)點(diǎn)的函數(shù)值。7、計算感興

4、趣的單元內(nèi)的函數(shù)值。5有限單元法簡介為什么要用有限單元方法（有限單元方法的優(yōu)勢）：1、有限單元方法可以對任意形狀的問題進(jìn)行靈活剖分，特別適 用于非常復(fù)雜的模型，最早應(yīng)用于結(jié)構(gòu)力學(xué)。2、有限單元方法是工程中應(yīng)用最廣泛的計算及數(shù)值模擬方法。3、所需要的復(fù)雜節(jié)點(diǎn)生成可以利用圖形界面軟件（如CAD）實(shí) 現(xiàn)。4、有很多商用有限單元軟件系統(tǒng)可以使用（如：ANSYS、 ADINA、SMART）。6有限單元法簡介有限單元方法的應(yīng)用：1、力學(xué)、航空、土木工程、汽車工程2、結(jié)構(gòu)分析（靜力學(xué)、動力學(xué)、線形、非線性）3、熱傳導(dǎo)和流體力學(xué)4、電磁學(xué)5、地質(zhì)力學(xué)6、生物學(xué)7有限單元法簡介有限單元方法在地球物理中的應(yīng)用：1

5、、地殼變形2、地球物理流體力學(xué) 地球動力學(xué) 地幔對流3、地電磁學(xué)4、波動傳播5、強(qiáng)震地面運(yùn)動、地震工程8有限單元法簡介有限單元方法在地球物理中的應(yīng)用：1、地殼變形2、地球物理流體力學(xué) 地球動力學(xué) 地幔對流3、地電磁學(xué)4、波動傳播5、強(qiáng)震地面運(yùn)動、地震工程9線性代數(shù)基礎(chǔ)線性方程組：其中x1、x2、xn是待求變量，以上方程組寫為矩陣形式：其中：10線性代數(shù)基礎(chǔ)列向量： 行向量：矩陣加減運(yùn)算： 其中：矩陣相乘運(yùn)算： 其中：其中：A：lxm矩陣，B：mxn矩陣，i=1,2,l，j=1,2,n一般地： 但是11線性代數(shù)基礎(chǔ)特殊矩陣：矩陣的轉(zhuǎn)置：對稱矩陣：單位矩陣： 并且有12線性代數(shù)基礎(chǔ)特殊的行列式：方

6、陣A的行列式是一個數(shù)，表示為 ，或 13線性代數(shù)基礎(chǔ)矩陣求逆：方陣A可逆的充分必要條件是它的行列式不為零，即 如果行列式為零，稱A為奇異矩陣。矩陣求逆：對一個方陣A，如果存在一個矩陣B，使得則B是A的一個逆矩陣。當(dāng)A可逆時，A的逆矩陣為：其中，A*為A的伴隨矩陣，由A的代數(shù)余子式組成，14線性代數(shù)基礎(chǔ)矩陣求逆：對逆矩陣有:線性方程組求解：對線性方程組 ，如果A可逆，則線性方程組求解的主要任務(wù)是求系數(shù)矩陣A的逆矩陣。常用求解技術(shù)，如：高斯消元法(Gaussian elimination)（對給定的線性方程組施行初等變換，將其變成一個同解的階梯形方程組，從而達(dá)到求解的目的。）迭代方法15線性代數(shù)基

7、礎(chǔ)正定矩陣：對于方陣A，如果對所有的非零的向量x，則A為正定矩陣。正定矩陣是非奇異矩陣。矩陣的微分、積分：設(shè)該矩陣對變量t的微分為：相應(yīng)的積分為：16有限單元法基礎(chǔ)有限單元方法最早是為了解決彈性靜力學(xué)問題而提出來的，因此一般按照方法的發(fā)展歷史介紹其基本概念，如單元(elements)、剛度矩陣(stiffness matrix)等等。需要的預(yù)備知識：胡克定律（Hookes Law）靜力平衡原理功的概念應(yīng)變能17有限單元法基礎(chǔ)一維問題：設(shè)有線性方程組設(shè)有向量y， 。一般地，方程組兩側(cè)同乘y不改變它的解： 考慮泊松方程：其中u是標(biāo)量場，f是源項(xiàng)，一維情況下：泊松方程兩側(cè)乘以任意函數(shù)v(x)， 18

8、有限單元法基礎(chǔ)一維問題：以上方程在整個求解區(qū)域D上求積分,為簡單起見，定義區(qū)域D為0，1，則積分為區(qū)域的離散化(discretization)：為了求出近似解，需要對區(qū)域進(jìn)行某種形式的離散化，有限單元法中，函數(shù)值只定義在離散化后的離散點(diǎn)上（這一點(diǎn)和有限差分法類似）：19有限單元法基礎(chǔ)一維問題：區(qū)域的離散化(discretization)：20有限單元法基礎(chǔ)一維問題：有限單元法的中心思想是用選取的基函數(shù)的線性組合來近似表示函數(shù)值：其中N為區(qū)域上的離散點(diǎn)數(shù)，ci為系數(shù)（實(shí)數(shù)）。如果基函數(shù) 選取合理，ci等于離散點(diǎn)i處的真實(shí)函數(shù)值，意味著在離散點(diǎn)i處的基函數(shù) ，其它點(diǎn)上的基函數(shù)值都為零。下面再看前面

9、的積分方程： 21有限單元法基礎(chǔ)一維問題：對該方程的左側(cè)作分部積分：設(shè)u在邊界處的微分為零，得到：原來的方程變?yōu)椋浩渲衭為待求的未知函數(shù)值。22有限單元法基礎(chǔ)一維問題：如果用近似值代替u，上式對近似值同樣成立：其中：是我們利用選取的基函數(shù)展開后給出的近似函數(shù)值。V是任意選取的實(shí)數(shù)值函數(shù)，如果上式對任意函數(shù)成立，它對下式顯然成立：因此對所有的基函數(shù)成立。23有限單元法基礎(chǔ)一維問題：將以上內(nèi)容綜合到一起：得到：24有限單元法基礎(chǔ)一維問題：系數(shù)ck為常數(shù)，所以對每一個基函數(shù)有：上式可以寫為矩陣形式：得到解為：bi是基函數(shù)的系數(shù)，對于特殊選取的基函數(shù)，這些系數(shù)值正好是離散點(diǎn)i上的函數(shù)值。25有限單元法

10、基礎(chǔ)一維問題：基函數(shù)：我們希望基函數(shù)具有如下性質(zhì)：可以選擇具有這種性質(zhì)的任何函數(shù)，最簡單的是如圖所示的線性函數(shù)（藍(lán)線：基函數(shù)值）26有限單元法基礎(chǔ)一維問題：基函數(shù)的梯度：為了構(gòu)造剛度矩陣(stiffness matrix)，我們需要計算基函數(shù)（藍(lán)線）的梯度（紅線）：27有限單元法基礎(chǔ)一維問題：剛度矩陣(stiffness matrix)：知道以上性質(zhì)的基函數(shù)后，我們可以計算矩陣形式方程中的Aij和gi：基函數(shù)是定義在區(qū)間0，1上的連續(xù)函數(shù)，28有限單元法基礎(chǔ)一維問題：剛度矩陣(stiffness matrix)：假設(shè)模型離散在等間隔節(jié)點(diǎn)上，設(shè):基函數(shù)可以寫為：29有限單元法基礎(chǔ)一維問題：剛度矩

11、陣(stiffness matrix)：基函數(shù)的梯度為：30有限單元法基礎(chǔ)一維問題：剛度矩陣(stiffness matrix)：方程中Aij的計算：31有限單元法基礎(chǔ)一維問題：剛度矩陣(stiffness matrix)：方程中Aij的計算：32有限單元法基礎(chǔ)一維問題：剛度矩陣(stiffness matrix)：方程中Aij的計算：33有限單元法基礎(chǔ)一維問題：邊界條件和源項(xiàng)：原始方程為：如果假設(shè)：最后變換為以下方程：加上邊界條件后，設(shè)：其中u(0)和u(1)是區(qū)域0，1邊界處的值。34有限單元法基礎(chǔ)一維問題：邊界條件和源項(xiàng)：原始方程為：加上邊界條件后，設(shè)：最后變換為以下方程：可以寫為矩陣形

12、式：35有限單元法基礎(chǔ)一維問題：邊界條件和源項(xiàng)：用圖形形象表示（系統(tǒng)通過修改后的源項(xiàng)來考慮邊界條件）：邊界條件源項(xiàng)邊界條件36有限單元法基礎(chǔ)一維問題：數(shù)值算例（等間距離散點(diǎn)）：偏微分方程：區(qū)域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);邊界條件：u(0)=u(1)=0Matlab 有限單元程序Matlab 有限差分程序37有限單元法基礎(chǔ)一維問題：數(shù)值算例（等間距離散點(diǎn)）：偏微分方程：區(qū)域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);邊界條件：u(0)=u(1)=0有限單元Matlab程序計算結(jié)果（藍(lán)色）有限差分Matlab程序計算結(jié)果（紅色）3

13、8有限單元法基礎(chǔ)一維問題：數(shù)值算例（等間距離散點(diǎn)）、邊界條件不為零：偏微分方程：區(qū)域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);邊界條件：u(0)=0.15 u(1)=0.05有限單元Matlab程序計算結(jié)果（藍(lán)色）有限差分Matlab程序計算結(jié)果（紅色）39有限單元法基礎(chǔ)一維問題：不等間距離散點(diǎn)時的剛度矩陣：40有限單元法基礎(chǔ)一維問題：數(shù)值算例（不等間距離散點(diǎn)）：偏微分方程：區(qū)域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);邊界條件：u(0)=u0 u(1)=u1剛度矩陣的計算程序41有限單元法基礎(chǔ)一維問題：數(shù)值算例（不等間距離散點(diǎn)）：偏微分

14、方程：區(qū)域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);邊界條件：u(0)=0.15 u(1)=0.0542有限單元法基礎(chǔ)總結(jié)：1、有限單元方法的分析中，我們用一組正交基函數(shù)的線性組合來給出定義在區(qū)域D上的函數(shù)的近似值，線性組合的系數(shù)相當(dāng)于某些離散節(jié)點(diǎn)上的函數(shù)值。2、對某些偏微分方程組，離散節(jié)點(diǎn)上的值可以通過求解線性方程組獲得，求解線性方程組得過程包括矩陣（有時是稀疏矩陣）的求逆。3、邊界條件在變換以后的方程中自然滿足，這是有限單元方法的優(yōu)勢之一（和有限差分法相比）43有限單元法-1D單元坐標(biāo)變換：我們希望用一組基函數(shù)的線性組合來近似給出區(qū)間a,b上的函數(shù)值u(x):其

15、中i是位置xi處的離散節(jié)點(diǎn)（單元的邊界）的編號。因?yàn)樗袉卧械幕瘮?shù)具有相同的形式，因此為簡單起見，我們將坐標(biāo)系統(tǒng)變換到一個局部的坐標(biāo)系統(tǒng)中：因此，單元定義在以下區(qū)間：44有限單元法-1D單元1D線性單元：對于1D直線單元，沒有形狀的選擇，但是單元的長度可以在區(qū)域內(nèi)變化。我們希望在每一個單元內(nèi)，函數(shù)u()用以下的線性函數(shù)來近似：離散節(jié)點(diǎn)定義在1,2=0,1處。我們要求：45有限單元法-1D單元1D線性基函數(shù)：我們利用離散節(jié)點(diǎn)1,2處的函數(shù)值給出了系數(shù)ci的值。接下來，我們利用離散節(jié)點(diǎn)上的值給出近似函數(shù)值：其中N1,2()是一維單元的線性基函數(shù)。46有限單元法-1D單元1D二次單元：我們希望在

16、每一個單元內(nèi)，函數(shù)u()用以下的二次函數(shù)來近似：離散節(jié)點(diǎn)定義在1,2,3=0,1/2,1處。我們需要：47有限單元法-1D單元1D二次基函數(shù)：同樣，我們可以用基函數(shù)的線性組合（系數(shù)為相應(yīng)的三個離散節(jié)點(diǎn)上的值）近似給出函數(shù)值：注意這里我們每一個單元里用到了三個離散節(jié)點(diǎn)。48有限單元法-1D單元1D三次基函數(shù)：用同樣的過程，我們可以推導(dǎo)出三次基函數(shù)的形式：49有限單元法-2D單元坐標(biāo)變換：現(xiàn)在我們討論二維（2D）問題的單元的幾何特征和基函數(shù)，同樣，我們希望在局部的坐標(biāo)系統(tǒng)中討論問題，以三角形為例：變換前變換后50有限單元法-2D單元坐標(biāo)變換：任何一個三個角位置為 （反時針順序）的三角形可以通過坐標(biāo)

17、變換變換為直角等腰三角形：如果=0，以上變化和1D坐標(biāo)變換相同。現(xiàn)在我們希望用線性形式給出函數(shù)的近似：和1D的推導(dǎo)過程類似，我們得到：51有限單元法-2D單元系數(shù)：系數(shù)可以通過線性方程組求解：矩陣A為：包含1D部分52有限單元法-2D單元線性基函數(shù)：根據(jù)矩陣A和線性方程組，可以計算出系數(shù)，給出三角形單元的線性基函數(shù)：53有限單元法-2D單元二次基函數(shù)：定義在三角形上的任意函數(shù)可以用二次函數(shù)近似為：在變換后的坐標(biāo)系里，形式為： 和1D類似，這時需要另外的三 個點(diǎn)，共有六個點(diǎn)54有限單元法-2D單元二次基函數(shù)：為了計算系數(shù)，我們計算每一個節(jié)點(diǎn)上的值： 矩陣A為：系數(shù)可以通過P點(diǎn)的值計算： 55有限

18、單元法-2D單元二次基函數(shù)：基函數(shù)為：56有限單元法-2D單元二次基函數(shù)：基函數(shù)為：57有限單元法-2D單元四邊形單元：坐標(biāo)變換：變換到局部坐標(biāo)系統(tǒng)中。58有限單元法-2D單元四邊形單元：線性基函數(shù)：得到矩陣A為：基函數(shù)為：59有限單元法-2D單元四邊形單元：二次基函數(shù)：得到矩陣A（8x8）和基函數(shù)，如：60有限單元法-2D單元總結(jié)：有限單元方法的基函數(shù)的求取規(guī)程：1、將坐標(biāo)系統(tǒng)轉(zhuǎn)換到局部的（單元的）坐標(biāo)系統(tǒng)；2、對在整個單元內(nèi)部取值的函數(shù)進(jìn)行線性（二次、三次）函數(shù)的近似；3、利用插值條件（特殊選取的基函數(shù)保證相應(yīng)點(diǎn)上的值為1，其它點(diǎn)上為0），得到系數(shù)值（離散節(jié)點(diǎn)上函數(shù)值的函數(shù)）4、利用求解線性方程組得到的系數(shù)推導(dǎo)出n個節(jié)點(diǎn)上的n個基函數(shù)。61有限單元法-例子聲波方程：如何利用有限單元法求解與事件有關(guān)的問題？聲波方程：其中v為聲波速度。用和上述類似的思想，方程兩側(cè)同乘以任意