《有限元素法》PPT課件.ppt

計算機在材料中的應(yīng)用,主講教師 武建軍,第四章 有限元素法,數(shù)值方法的本質(zhì)就是把計算域內(nèi)有限位置（節(jié)點）上的因變量值當(dāng)作基本未知量來處理，任務(wù)是提供一組關(guān)于這些未知量的代數(shù)方程并規(guī)定求解這些方程的算法。 有限元法是以變分原理為基礎(chǔ)吸取差分格式的思想而發(fā)展起來的一種有效的數(shù)值解法。 有限元素法（FEM）是求解微分方程的一個常用數(shù)值解法。,FEM的應(yīng)用領(lǐng)域,50年代，有限元法首先從結(jié)構(gòu)力學(xué)中得到應(yīng)用，而后推廣到流體力學(xué)、熱傳導(dǎo)等領(lǐng)域,現(xiàn)在已被廣泛用于物理和工程設(shè)計計算的很多領(lǐng)域。該方法逐漸成為求解許多數(shù)學(xué)物理方程邊值問題的高效能的方法。,1 FEM的基本思想及其特點,利用FEM方法求解問題包括單元劃分、假定單元分布函數(shù)、列出單元方程、聯(lián)立求解等幾個步驟。 例如解決一個超靜定梁的問題，每一個桿件就可以作為一個單元來處理,可以列出自己的載荷與位移的關(guān)系方程(單元方程),然后聯(lián)立求解。,1-2 FEM的基本思想-里茲法,用插值（或其它）方法用節(jié)點處的函數(shù)值構(gòu)建單元分布函數(shù)； 將代入原來的微分方程，將得到有關(guān)節(jié)點變量值的代數(shù)方程組。 求解代數(shù)方程組，獲得節(jié)點函數(shù)值。 離散化方程仍受未知函數(shù)的微分方程支配。,單元劃分的意義,在整個計算域內(nèi)尋找試探函數(shù)不太容易。如果將計算域劃分成一定數(shù)量的單元，單元分布函數(shù)容易建立。 因為選取分段分布的結(jié)果，一定的離散方程只與少數(shù)幾個節(jié)點有關(guān)。一方面分布函數(shù)容易構(gòu)造，而且也容易求解。,分布函數(shù)的隨意性,由于分布假設(shè)以及推導(dǎo)方法的不同，一個微分方程的離散化方程不是惟一的。 因為只需要滿足本質(zhì)性邊界條件，而不必考慮自然邊界條件（第二、第三類邊界條件自動滿足），試探函數(shù)的選取是比較容易的。 試探函數(shù)階次提高，解的精度也提高。,當(dāng)網(wǎng)格特別細(xì)密時，相鄰節(jié)點之間的變化就很小，因此單元內(nèi)分布假設(shè)的實際細(xì)節(jié)變得不再重要。離散化方程的解將趨近于相應(yīng)微分方程的精確解。,單元形狀,單元的形狀沒有限制。例如橢圓孔應(yīng)力集中集中問題，圖中將它劃分為三角形網(wǎng)絡(luò)。把原來的連續(xù)體簡化為由有限個三角形單元組成的離散體。其中三角形單元之間只在節(jié)點處用鉸鏈相連，把載荷按照靜力等效原則也轉(zhuǎn)移到節(jié)點處。,1-2 FEM的特點,物理概念清晰，特別是對于力學(xué)問題。 靈活性與通用性。由于單元形狀靈活，易于處理復(fù)雜區(qū)域、復(fù)雜邊界條件。 而對于具有規(guī)則的幾何特性和均勻的材料特性問題,差分法的程序設(shè)計比較簡單,收斂性也比有限元法好。 有限元法同時具有里茲法與差分法的優(yōu)點，使變分問題的直接解法變成了工程計算中的現(xiàn)實。,FEM的特點,有限元素方法是物理量的矩陣分析方法在連續(xù)體中的有效推廣。每個元素都采用有限個參數(shù)來描述它的物理特性。 有限元素方法是基于虛功原理,或者說是變分原理。它不象差分法那樣直接去解場方程,而是求解一個虛功取極小值的變分問題。 FEM是解決復(fù)雜區(qū)域、邊界條件數(shù)學(xué)物理方程邊值問題的一種比較完美的離散化方法。,本章提要,有限元素法的原理概要 舉例說明如何運用有限元素法 彈性力學(xué)平面問題 熱傳導(dǎo)問題 比較有限元素法和有限差分法。,2 彈性力學(xué)基礎(chǔ)知識,彈性力學(xué)研究宏觀均勻、各向同性固體的彈性變形，例如刃型位錯應(yīng)力場計算。 嚴(yán)格地說，任何彈性體總是處在空間應(yīng)力狀態(tài)，因而實際問題都是三維空間問題。但是，有些彈性力學(xué)問題可以簡化為平面問題。,平面應(yīng)力問題,例如平面薄板的拉伸變形問題，由于厚度很小，而載荷又平行板面且沿厚度方向均勻分布，因此可以近似認(rèn)為沿厚度方向的應(yīng)力分量等于零。,平面應(yīng)變問題,水庫大壩的長度比高度和寬度要大得多，而載荷又都與橫斷面平行且沿長度方向均勻分布，可以認(rèn)為沿長度方向的應(yīng)變分量等于零，這種問題稱平面應(yīng)變問題。,2-1 變形參數(shù),單元體彈性變形參數(shù)包括 位移沿x,y,z軸的三個分量 u,v,w 力或載荷沿x,y,z軸的三個分量 X,Y,Z 應(yīng)變ij（作用面垂直于i軸，指向j方向）i,j=x,y,z 應(yīng)力ij(規(guī)定同應(yīng)變 ),一般規(guī)定,當(dāng)i=j, ij為正應(yīng)變,表示線段伸長或縮短,可簡化為i ，如x。一般規(guī)定，伸長應(yīng)變?yōu)檎怠?ij， ij表示切應(yīng)變（角應(yīng)變），表示兩線段之間夾角的變化。一般規(guī)定，直角變成銳角切應(yīng)變大于零。 應(yīng)變分量有9個,一般有6個獨立分量。 應(yīng)力分量與應(yīng)變分量類似，也有9個。,2-2 幾何方程,2-3 廣義虎克定律,廣義虎克定律,拉梅方程的其它形式,廣義虎克定律,平面應(yīng)力問題 的虎克定律,平面應(yīng)變問題的虎克定律,2-4 力學(xué)平衡方程,應(yīng)力與體積力，如重力之間的平衡 其中fi分別是體積力在i方向上的分量。,力學(xué)平衡方程,應(yīng)力與表面力之間的平衡 式中X、Y、Z是表面力的三個分量，l,m,n是表面外法線的方向余弦。,3 變分方法與虛功原理,有限元法建立離散方程的方法有三類。 直接法 例如，超靜定桁架問題，每個組件就是一個元素。易于理解,但只能用于較簡單的問題，實際用途不大。 變分法 把有限元法歸結(jié)為求泛函的極值問題。使有限元法建立在更加堅實的數(shù)學(xué)基礎(chǔ)上,擴大了有限元法的應(yīng)用范圍。 加權(quán)余量法 直接從基本微分方程出發(fā),求得近似解。對于根本不存在泛函的工程領(lǐng)域都可采用,從而擴大了有限元法的應(yīng)用范圍。,3-1 變分原理,變量與變量間的關(guān)系稱為函數(shù)。 如果對于某一類函數(shù)y(x)中的每一個函數(shù)y(x),I都有一值與其對應(yīng)，則變量I叫做依賴于函數(shù)y(x)的泛函記為Iy(x)。 泛函是函數(shù)集合的函數(shù),也就是函數(shù)的函數(shù)。,泛函舉例-曲線長度,在直角坐標(biāo)系中兩定點M和N,連接兩點曲線的長度L與曲線的形狀y=y(x)有關(guān)。曲線長度是由曲線的方程y=y(x)所確定,因此，曲線的長度L就是曲線函數(shù)的一個泛函，可以記為Ly(x)。,泛函舉例-應(yīng)變能,彈性體在各種力作用下引起彈性應(yīng)變（x）, 物體具有不同的位移函數(shù)（或應(yīng)變分布函數(shù)）具有不同的應(yīng)變能，那么應(yīng)變能就是一個依賴位移函數(shù)的泛函。,其它例子,外力勢能F也是位移函數(shù)的泛函 f為體積力，q為面積力。 此外總勢能P=V-F也是位移函數(shù)的泛函。,變分的定義,如果對應(yīng)于函數(shù)y(x)的微小改變，有泛函的微小改變與之對應(yīng)，稱泛函是連續(xù)的。 函數(shù)的變分y是指兩個非常接近的函數(shù)y(x)與y1(x)的差 y= y(x)y1(x) 而泛函的變分是,變分問題,變分問題就是研究泛函的極值問題。這個問題的解法類似于函數(shù)極值的求法。 函數(shù)f(x)的極值可以根據(jù)df=0這個條件去尋找判斷。類似于普通函數(shù)取極值的條件,若具有變分的泛函在y=y0(x)處取得極值,那末泛函在該處的變分應(yīng)為零。用數(shù)學(xué)公式寫出為 即可以根據(jù)變分等于零這個條件去尋找y0(x)。,一維歐拉方程,可以證明，泛函Q 實現(xiàn)極值的必要條件是函數(shù)y(x)滿足一維歐拉方程,二維歐拉方程,若f(x,y)在R內(nèi)二階可微，區(qū)域R的邊界分為B、C兩部分，在B邊界滿足f=f1(x,y)，C邊界上為未知函數(shù)。泛函 函數(shù)G（f）沿邊界C取值。 假定泛函有極值，必須滿足歐拉方程 在區(qū)R內(nèi)， 在邊界上， 式中l(wèi),m是邊界外法線的方向余弦。,變分問題舉例,變分原理,由于使泛函實現(xiàn)極值的函數(shù)必然滿足相應(yīng)的歐拉方程，這樣可以把求解某些微分方程的問題轉(zhuǎn)化為泛函極值問題。如果求得了使泛函實現(xiàn)極值的函數(shù)，也就是微分方程的解函數(shù)。 變分方法的主要難點就是尋找適當(dāng)形式的泛函（許多微分方程不存在對應(yīng)的泛函）。泛函的尋找一般要根據(jù)物理概念，如虛功原理、最小能量原理等來進行。,彈性力學(xué)問題的變分原理,最小總勢能原理指出，在滿足邊界條件的所有位移中，以保持力的平衡狀態(tài)的位移造成的總勢能最小。 因此，求平衡態(tài)的位移函數(shù)就等同于求應(yīng)變能（位移函數(shù)的泛函）的極小值，即P(u0(x)=0的解。這就是彈性力學(xué)問題有限元方法的變分原理。,當(dāng)單元內(nèi)任意一點的函數(shù)值根據(jù)插值法用節(jié)點函數(shù)值表示時，單元積分也成為節(jié)點函數(shù)值的函數(shù)。 由于求解域R已經(jīng)劃分成若干單元，泛函變成這些子域上積分的和。顯然，泛函也成為節(jié)點函數(shù)值 （i=1,2,3n）的函數(shù)。 由泛函極值的必要條件，得方程組 ，i=1,2,3,n,3-2 虛功原理,設(shè)變形體處于平衡受力狀態(tài)，體積力為f,在自由邊界上的表面力為q。設(shè)變形體產(chǎn)生任一虛位移u*，相應(yīng)的虛應(yīng)變?yōu)?=Bu*，則體積力和表面力所作的虛功，恒等于應(yīng)力在虛應(yīng)變上所作的虛變形功（應(yīng)變能增量）,如果外力是集中力P1,P2，P3Pn,而各力相應(yīng)的虛位移分量*i, 以上結(jié)果與變分原理的結(jié)果相同。,3-3 加權(quán)余量法,用變分法求解微分方程，首先要找到相應(yīng)的泛函。 對于有些問題相應(yīng)的泛函尚未找到，或者根本不存在相應(yīng)的泛函。在這種情況下，就無法用變分法求解。 加權(quán)余量法（也稱加權(quán)余值法）是求微分方程近似解的一種有效方法。,設(shè)有微分方程 假設(shè)有一個滿足邊界條件和具有一定連續(xù)程度的試探函數(shù) （其中含有若干待定系數(shù)）使 余量R與所選的試探函數(shù)有關(guān)，希望R在某種意義上比較小。通過把余量與加權(quán)函數(shù)Wi(x)正交化的途徑，化為代數(shù)方程組，即令,加權(quán)函數(shù),比較常用的加權(quán)函數(shù)有 狄拉克函數(shù)(x),配置法； 冪函數(shù) ,i=0,1,2,矩量法； 里茲基函數(shù)，Galerkin法；精度高，常用 余數(shù)R,最小二乘法。 有限元法中積分在單元內(nèi)進行。,4 彈性力學(xué)問題的FEM,單元分析 總體分析 邊界條件處理 方程求解,4-1 單元分析,單元分析的任務(wù)是推導(dǎo)基本未知量節(jié)點位移與節(jié)點力F之間的關(guān)系。 一般可從節(jié)點位移、單元內(nèi)各點位移(插值)、應(yīng)變（幾何方程）、應(yīng)力（虎克定律）入手，最終求得節(jié)點位移與節(jié)點力F之間的關(guān)系,平面三角形單元,圖示一個三角形單元 三個節(jié)點按逆時針方向的順序編碼為i,j,m。 節(jié)點坐標(biāo)分別為（xk,yk）k=i,j,m。,在彈性力學(xué)平面問題中，每個節(jié)點有兩個位移分量，因此三角形單元的位移向量可寫成 與此對應(yīng)的物理量是六個節(jié)點力分量 單元分析的任務(wù)就是推導(dǎo)以上兩向量之間的關(guān)系，即 其中轉(zhuǎn)換矩陣Ke稱單元剛度矩陣。,位移插值函數(shù),單元（元素）分析的第一步是由單元節(jié)點位移分量推算單元內(nèi)部任意一點的位移。在任一單元內(nèi)可以通過線性插值構(gòu)造三角形內(nèi)部任意一點的位移u和v，即可設(shè)在,代入節(jié)點位移可得到方程組，解方程組得到系數(shù)值,其中 按ijmi的順序輪換可得到其它參數(shù)。三組常數(shù)都是與節(jié)點坐標(biāo)有關(guān)的常數(shù). A是三角形元素的面積 ak,bk,ck是行列式第k行對應(yīng)元素的代數(shù)余子式。,位移插值函數(shù),將系數(shù)i代入u,v表達式，整理可得 式中N是插值基函數(shù),也稱形函數(shù) l=i,j,m。,形函數(shù)的性質(zhì),形函數(shù)描述了單元內(nèi)部研究對象例如位移函數(shù)的分布情況。 在單元任意一點，三個形函數(shù)的和等于一（設(shè)三節(jié)點位移相等，就可以證明）。 在單元的邊上（如ij邊），有一個形函數(shù)（如Nm）等于零。 在單元節(jié)點上,單元剛度矩陣,（1）由位移求應(yīng)變 幾何方程用矩陣表示，B稱幾何矩陣。由于插值函數(shù)是一次函數(shù)，在單元內(nèi)的應(yīng)變分量都是常量。,（2）位移與應(yīng)力 根據(jù)虎克定律， 其中轉(zhuǎn)換矩陣S=DB (3)由節(jié)點位移求節(jié)點力 根據(jù)虛功原理 ，,單元剛度矩陣的性質(zhì),1 單元剛度矩陣是66對稱矩陣,符號規(guī)定,剛度矩陣中kpq,p代表節(jié)點力編號和方向， q代表位移的編號與方向 帶代表水平方向，不帶代表垂直方向。 例如， 代表m節(jié)點力垂直方向分量Vm與i節(jié)點水平位移分量ui之間的作用系數(shù)。,2 K是奇異矩陣，對應(yīng)的行列式等于零，不存在逆矩陣。實際無法獨立求解 從幾何關(guān)系也可得到這一結(jié)論。如果給定節(jié)點位移，按照這個關(guān)系可以計算惟一的節(jié)點力。但是反過來，根據(jù)節(jié)點力卻不能得出節(jié)點位移的唯一解。因為單元可以產(chǎn)生任意的剛體位移-不影響節(jié)點力。,4-2 整體分析,平面彈性力學(xué)有限元法計算的最終結(jié)果是要求出區(qū)域R中的位移分布、應(yīng)變分布及應(yīng)力分布。 由于已經(jīng)把解域R劃分成有限個三角形單元及若干個節(jié)點，并把整個場離散到這些節(jié)點上，因此，我們的任務(wù)就是把這些節(jié)點上的位移、應(yīng)變和應(yīng)力求出來。,整體分析的任務(wù),由于一個節(jié)點不止屬于一個單元，因而要將單元方程加以合并組合，得到反映所有節(jié)點節(jié)點力與位移關(guān)系的總方程組。即根據(jù)單元剛度矩陣建立連續(xù)體的總體剛度矩陣、根據(jù)邊界條件建立節(jié)點載荷向量。 整體分析的任務(wù)就是建立整體剛度矩陣和節(jié)點載荷向量。,合成方法與步驟,在剛度矩陣合成時，存在兩個問題。一是單元矩陣與總體矩陣的階數(shù)不同；二是單元節(jié)點是局部號碼。 剛度矩陣的集成分兩步進行。第一步把單元剛度矩陣擴大到nn階，然后把單元剛度矩陣中Ke的9個2 2子塊搬家，變成按總體編碼的9個子塊,其它子塊用零來填充。例如某單元局部編碼i,j,m對應(yīng)于總碼（節(jié)點編號）2,4,5,則該單元矩陣在擴大、搬家之后將變成（設(shè)n=7）,擴維,疊加,第二步，把擴充以后的單元矩陣迭加，即得到整體剛度矩陣。為了節(jié)省存儲容量，以上兩步實際上是穿插進行的，即按對號入座、邊搬家、邊累加的方法進行集成。 節(jié)點載荷向量的構(gòu)建方法與剛度矩陣相似。由作用力與反作用力的關(guān)系，大多數(shù)節(jié)點的節(jié)點力的合力都等于零，只有那些第一類邊界節(jié)點或受到體積力的節(jié)點才不等于零。,4-3 第一類邊界條件的處理,在有限元法中，自然邊界條件（包括第二 、第三類邊界條件）已經(jīng)包含在變分公式當(dāng)中，因而是自動滿足的。 第一類邊界條件必須強制性地引入到有限元方程當(dāng)中，即對有關(guān)邊界節(jié)點的方程進行修改。為了修改方便，總體節(jié)點編號一般先編內(nèi)部節(jié)點，最后編第一類邊界節(jié)點。,例如節(jié)點9的水平位移 ,只需要將K中第17行主對角元素改為1,其余全為零，對應(yīng)的載荷向量也改成零。 例如節(jié)點7的垂直位移 ,則可進行如下修改。將K中第14行主對角元素乘以一個較大的數(shù)M（例如108），對應(yīng)的載荷向量也改成 ,其余各項不變。這時，這個方程中除包含大數(shù)M的兩項外，其余各項相對較小，可以忽略不計，可以得到 。,4-4 方程特點與求解,1.整體剛度矩陣的特點 （1）整體剛度矩陣是一個對稱矩陣。 （2）整體剛度矩陣是一個稀疏矩陣，它的大多數(shù)元素都是零。 （3）整體剛度矩陣中的非零元素分布在以主對角線為中心的斜帶型區(qū)域內(nèi)。 利用矩陣的對稱性和帶型分布規(guī)律，在計算機中可以只存儲上半帶元素。,2.方程的求解 可以用帶消去法或高斯賽德爾迭代法等求解方程組，得到節(jié)點位移。 然后再用單元分析公式計算每個單元的應(yīng)變或應(yīng)力。,5 彈性力學(xué)問題有限元解法,下面以一個簡單的模型，示范彈性力學(xué)平面問題有限元解法的基本程序。 如圖所示的試件，受到一組集中力。,5-1 單元分析,第一單元，節(jié)點坐標(biāo)分別是i=1(0,20),j=2(0,10),m=3(10,10)（逆時針方向編號），面積A=50,彈性系數(shù),設(shè)為平面應(yīng)力問題，泊松比0.28, E=21000，板厚t=10,1，2，4單元剛度矩陣,對于第二單元，節(jié)點坐標(biāo)分別是i=2(0,10),j=4(0,0),m=5(10,0),A=50 B及K與第一、四單元相同,第三單元，節(jié)點坐標(biāo)分別是 i=2(0,10), j=5(10,0), m=3(10,10) A=50,單元3剛度矩陣,5-2 整體分析,第三單元剛度矩陣擴容后的形式,整體剛度矩陣,強加第一類邊界條件后剛度矩陣,5-3 平衡方程及其解,其中，,節(jié)點位移,解方程，可以得到節(jié)點位移。,在單元內(nèi) 利用位移與應(yīng)變的關(guān)系（幾何方程）可得到節(jié)點應(yīng)變； 根據(jù)虎克定律可以計算節(jié)點應(yīng)力。 例如對于第一單元，根據(jù)節(jié)點位移，計算1 2 3三節(jié)點的應(yīng)變和應(yīng)力,節(jié)點應(yīng)變,節(jié)點應(yīng)力,6二維熱傳導(dǎo)問題FEM解法,插值函數(shù)中的系數(shù),A為三角形單元的面積,因為傳熱問題沒有對應(yīng)的泛函，以下我們采用加權(quán)余量法構(gòu)造單元方程。 應(yīng)用加遼金法，令,因為所取的插值函數(shù)為一次函數(shù)，首先應(yīng)用分部積分公式，把二階導(dǎo)數(shù)項化為一階導(dǎo)數(shù) S是單元e的周界。,如果單元周界外法線n與x軸的夾角余弦為nx, 與y軸的夾角余弦為ny,上式中沿單元周界的積分項可以寫成，,此時，微分方程變成如下形式 如果溫度插值函數(shù)使得溫度在每個單元之間連續(xù)的話，則經(jīng)過單元集合，相鄰邊界的積分必將互相抵消，上式中的單元積分僅在邊界單元才有作用。,單元溫度分布,單元方程,對于內(nèi)部單元，不必計算曲線積分 對于熱流邊界單元,換熱邊界單元方程,規(guī)定記號,6-2 單元分析,對于常用的三角形線性單元，由于 l=i,j,m 系數(shù)都是節(jié)點坐標(biāo)的函數(shù)。一階導(dǎo)數(shù) 均為常數(shù)，可以提到積分號之前,單元熱傳導(dǎo)矩陣,單元熱交換矩陣,根據(jù)里茲基函數(shù)的性質(zhì),單元熱容矩陣,向量,單元方程,對于內(nèi)部單元 對于熱流邊界單元 對于換熱邊界,整體傳熱方程,把所有單元集合可求得總體方程 其中,6-3 時間域的離散,不穩(wěn)定溫度場計算含有 是一個常系數(shù)微分方程，可以用差分方法求得數(shù)值解。,向后差分格式,將 代入 得到,C-N格式,加遼金格式,兩點格式通式,三點向后差分公式,式中a是變步長系數(shù)。如果a=1，步長固定。 經(jīng)過時間域的離散化，最后可將問題歸結(jié)為求解代數(shù)方程組問題。,7 FEM與FDM的比較-處理方法,有限元素法是基于變分原理,同時吸收了差分法中的區(qū)域劃分的思想。它將所要求解的物理問題化為泛函求極值的問題,再通