有限元法基礎

1、有限元法基礎華中科技大學材料成形與模具技術國家重點實驗室教材n董湘懷. 材料成形計算機模擬，第2版. 北京：機械工業(yè)出版社，2005 n講課：24學時主要參考書n王勖成. 有限單元法. 北京：清華大學出版社，2003 nT. R. Chandrupatla, A.D. Belegundu著，曾攀譯. 工程中的有限元方法. 北京：清華大學出版社，2006n曾攀. 有限元分析及應用. 北京：清華大學出版社，2004主要參考書n謝水生，李雷. 金屬塑性成形的有限元模擬技術及應用. 北京：科學出版社，2008n李尚健. 金屬塑性成形過程模擬. 北京：機械工業(yè)出版社，1999 n鐘志華，李光耀. 薄板沖

2、壓成型過程的計算機仿真與應用. 北京：北京理工大學出版社，1998材料成形問題的復雜性n非線性 幾何非線性 物理非線性 狀態(tài)非線性n模具和工藝的復雜性 復雜幾何形狀 多工序n多物理場耦合 變形 熱傳導n材料組織性能變化 相變、再結晶 織構 損傷數(shù)值模擬方法的優(yōu)越性n經(jīng)濟、快速、優(yōu)化、并行n 結果詳盡 應力、應變、溫度 組織性能變化n虛擬、靈活板料成形數(shù)值模擬實例板料成形數(shù)值模擬實例有限元法概述n有限元法的發(fā)展過程 1960年，Clough求解平面彈性問題，第一次提出“有限單元法”的名稱n有限元法的工程應用 桿件結構問題 平板彎曲問題 殼體問題 熱傳導問題 動力學問題 流固耦合問題 接觸碰撞問題

3、有限元法概述n有限元分析軟件 NASTRAN，SAP， MARC，ANSYS，ABAQUS， JIFEX， LS-DYNA3D n有限元法的發(fā)展方向 本構模型，單元模型，數(shù)值分析方案，CAD/CAE/CAM集成 有限元法的基本思想n 化整為零n 積零為整有限元法的基本思想n基本思想 1）將連續(xù)的求解系統(tǒng)離散為一組由節(jié)點相互聯(lián)在一起的單元組合體 2）在每個單元內(nèi)假設近似函數(shù)來分片表示系統(tǒng)的求解場函數(shù) 有限元法的基本思想n有限元法分類 1）位移法：基于最小勢能原理或虛功原理 2）力法： 基于最小余能原理 3）雜交法：基于修正余能原理 4）混合法：基于Reissner變分原理 有限元法的基本思想n位

4、移法基本過程 1）離散化過程 3）約束處理過程 2）單元平衡方程組裝過程 5）應變、應力回代過程 4）方程組求解過程 離散化過程n最小勢能原理 P V A G 彈性體 彈性體的勢能peipWW 為彈性體變形后所具有的內(nèi)能 iW為彈性體所受的外力功 eWdVVT21VAdVdAGuPuTT離散化過程 為彈性體的應變 為彈性體的應力 u為彈性體的可容位移 彈性體處于平衡狀態(tài)時，其勢能應為最小 P0TTTVAVdVdAdVGuPu0離散化過程n單元插值關系 eNuu n單元幾何關系 n單元本構關系 Lu DeN為單元形函數(shù)矩陣 L為單元幾何微分算子 為單元彈性矩陣 eD0TTTvavPdvdadvG

5、uPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDBeu為單元節(jié)點自由度向量離散化過程0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDB0TTTvavPdvdadvGuPuB 稱為應變矩陣 LNB fku e單元平衡方程或單元剛度方程 k 稱為單元剛度矩陣 vedvBDBkT f 稱為單元載荷向量 vadvdaGNPNfTT單元剛度矩陣的特性 n對稱性 n奇異性 n主元恒正且對角占優(yōu) 離散化過程線彈性問題幾何方程三維問題 Luwvuxzyzxyzy

6、xzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000三維問題Lu Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000線彈性問題幾何方程二維問題 Luvuxyyxyvxuyvxuxyyyxx00二維問題平面應力和平面應變狀態(tài) 線彈性問題幾何方程二維問題 二維問題軸對稱狀態(tài) Luwurzzrrrwzuzwrururzzzrr0010Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000線彈性問題幾何方程一維問題 一維問題 LuuxxuxxLuvuxyyxyvxuy

7、vxuxyyyxx00線彈性問題本構方程三維問題 三維問題De221000000221000000221000000100010001)21)(1 (EeDE為彈性模量；為泊松比 線彈性問題本構方程平面應力 二維問題平面應力狀態(tài) 00yzxz000yzxzzzzxyzxyzzyyxx000 xyyyxxxyyyxxzxyzxyzzyyxxxyzzyyxx線彈性問題本構方程平面應力 zxyzxyzzyyxxezxyzxyzzyyxxD2100010112EeDxyyyxxxyyyxxE2100010112平面應力狀態(tài) 00000 xyzzyyxxexyyyxxD線彈性問題本構方程平面應變 二維問

8、題平面應變狀態(tài) 000yzxzzz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyyyxx00000 xyyyxxexyzzyyxxD線彈性問題本構方程平面應變 zxyzxyzzyyxxezxyzxyzzyyxxDxyyyxxxyyyxxE221000101)21)(1 (平面應變狀態(tài) 221000101)21)(1 (EeD線彈性問題本構方程軸對稱 二維問題軸對稱狀態(tài) 00yzxz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyzzyyxx線彈性問題本構方程軸對稱 二維問題軸對稱狀態(tài)

9、00zr00zrzrzrzzrrzrzzrr00zrzzrrzrzrzzrrzrzzrr線彈性問題本構方程軸對稱 軸對稱狀態(tài) zxyzxyzzyyxxezxyzxyzzyyxxDzrzrzzrrezrzrzzrrDzrzzrrezrzzrr0000D221000010101)21)(1 (EeDzrzzrrzrzzrrE221000010101)21)(1 (線彈性問題本構方程一維問題 一維問題xxxxEEeD常用單元模型 n單元模型插值關系一一對應n單元類型一維單元、二維單元、三維單元等參單元、超參單元、次參單元常用單元模型n一維單元 2節(jié)點線單元 1 2 1 2 3 1 2 3節(jié)點線單元梁

10、單元常用單元模型n二維單元3節(jié)點三角形線性單元 1 2 3 6節(jié)點三角形二次單元 1 2 3 5 6 4 常用單元模型n二維單元10節(jié)點三角形三次單元 10 1 2 3 6 9 4 5 7 8 4節(jié)點四邊形雙線性單元 1 2 3 4 常用單元模型n二維單元8節(jié)點四邊形二次單元12節(jié)點四邊形三次單元 1 2 3 4 8 7 6 5 8 7 1 2 3 4 11 9 5 6 10 12 常用單元模型n三維單元4節(jié)點四面體線性單元10節(jié)點四面體二次單元 1 2 3 4 1 10 9 8 4 7 2 3 6 5 常用單元模型n三維單元8節(jié)點六面體線性單元20節(jié)點六面體二次單元 8 6 1 2 3 4

11、5 7 8 16 17 10 3 20 19 18 15 14 6 13 12 11 9 2 4 1 5 7 常用單元模型n準三維空間單元桁架單元一維2節(jié)點線單元+單元局部隨體坐標系 為什么要建立單元局部隨體坐標系 ？1.簡化分析問題的復雜程度。2.在局部坐標系中，空間桁架的每根桿都變成了一維2節(jié)點線單元常用單元模型n準三維空間單元框架單元三維梁單元+一維2節(jié)點線單元+單元局部隨體坐標系 兩端都是剛性聯(lián)結 可以承受拉壓、彎曲、扭轉(zhuǎn)3種變形模式 框架單元的特點常用單元模型n準三維空間單元板單元薄板單元中厚板單元彎曲和橫向剪切2種變形模式抵抗板的變形如果板很薄，忽略橫向剪切抗力，認為抵抗載荷的主要

12、因素是彎矩常用單元模型n準三維空間單元殼單元 抵抗拉壓變形的二維單元+板單元+單元局部隨體坐標系。適合于薄殼單元和中厚殼單元從幾何上分為薄殼單元和中厚殼單元 組合單元常用單元模型n準三維空間單元 殼理論單元 由空間殼理論嚴格構造的殼單元。適合于薄殼單元和中厚殼單元 退化單元 由三維實體單元退化成的殼單元。只適合于中厚殼單元 單元模型構造 n有限元法的基本思想 通過單元分片近似，在每個單元內(nèi)假設近似函數(shù)來分片表示系統(tǒng)的場函數(shù) n選擇近似函數(shù)簡單、實用的原則在有限元法中，近似函數(shù)稱為插值函數(shù) 單元模型構造n插值函數(shù) 一般都采用多項式函數(shù)，主要原因是: 采用多項式插值函數(shù)比較容易推導單元平衡方程，特

13、別是易于進行微分和積分運算。隨著多項式函數(shù)階次的增加，可以提高有限元法的計算精度。從理論上說，無限提高多項式的階數(shù)，可以求得系統(tǒng)的精確解。單元模型構造方法 n整體坐標系法n局部坐標系法 nLagrange插值方法nHermite插值方法n1）Lagrange插值方法 單元模型構造基本思想與方法 n對于有n個節(jié)點的單元，如果它的節(jié)點參數(shù)中只含有場函數(shù)的節(jié)點值ui，i=1,2,n ，則單元內(nèi)某一確定自由度方向的場函數(shù)u可插值為 1niiiuN un式中，Ni是插值函數(shù)，它有以下性質(zhì) 單元模型構造基本思想與方法 1(,)1ijjjijniiN xyzNn式中，ij是kronecker符號，xj、yj

14、、zj 是單元節(jié)點坐標值。 n上式也反映了插值函數(shù)的一個共同特性：插值函數(shù)Ni在i節(jié)點處的值等于1，在其他節(jié)點處的值于0，也稱為插值函數(shù)的正交性。 單元模型構造基本思想與方法 n插值函數(shù)的和等于1，也稱為插值函數(shù)的正規(guī)性。 n2）Hermite插值方法 單元模型構造基本思想與方法 n如果希望在單元間的公共節(jié)點上還保持場函數(shù)導數(shù)的連續(xù)性，則節(jié)點參數(shù)中還應包含場函數(shù)導數(shù)的節(jié)點值。這時可以采用Hermite多項式作為單元的插值函數(shù)。例如，對于只有兩個節(jié)點的一維單元，形函數(shù)Ni采用Hermite多項式插值為 單元模型構造基本思想與方法 221(0)11d( )( )( )diiiiiiuuHuH或者4

15、1( )( )iiiuNQ式中(0)(0)1122(1)(1)3142( ),( )( ),( )NHNHNHNH11223412,dd,ddQu QuuuQQnHermite多項式具有以下性質(zhì) 單元模型構造基本思想與方法 (0)(0)(1)(1)d( )(), 0dd( )()0, djjiijijiijijHHHHn上式在兩個節(jié)點處最高保持場函數(shù)的一階導數(shù)連續(xù)性，這種多項式稱為一階Hermite多項式。零階Hermite多項式就退化為Lagrange多項式。進一步推廣，在節(jié)點處保持至場函數(shù)的n階導數(shù)連續(xù)性，就稱為n階Hermite多項式。 單元模型構造基本思想與方法 221(0)11d(

16、)( )( )diiiiiiuuHuH單元模型構造方法n2節(jié)點線單元12 oxu1u2x1x2ux1. 假設插值多項式xaaxu10)(2. 利用節(jié)點值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua單元模型構造方法3. 代入a0 和 a1，得插值多項式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同類項，設 l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu單元模型構造方法n關鍵 如何構造插值多項式 u ?二維問題三維問題，如何構造插值多項式?n收斂性條件

17、在單元內(nèi)，場函數(shù)必須是連續(xù)的； 完備性：插值多項式的階次必須由低到高依次增加，不能出現(xiàn)跳躍現(xiàn)象； 協(xié)調(diào)性：各單元邊界必須連續(xù)，單元邊界不能出現(xiàn)開裂現(xiàn)象。 插值多項式收斂性條件 n收斂：當單元逐漸縮小時，如果插值多項式滿足收斂性條件，則數(shù)值解將收斂于精確解 插值多項式收斂性條件n協(xié)調(diào)單元 滿足插值多項式收斂性條件和的單元 n完備單元 滿足插值多項式收斂性條件的單元ncr 階連續(xù)性 插值多項式的第r階導數(shù)是連續(xù)的 插值多項式收斂性條件n非協(xié)調(diào)單元與部分協(xié)調(diào)單元 對于一般固體力學問題來說，協(xié)調(diào)性要求單元在變形時，相鄰單元之間不應引起開裂、重疊或其它不連續(xù)現(xiàn)象。例如，梁、板、殼等單元，在單元邊界不但要

18、求位移是連續(xù)的，而且其一階導數(shù)也必須是連續(xù)的。板、殼單元位移函數(shù)沿單元邊界的法向?qū)?shù)（轉(zhuǎn)角）的連續(xù)性一般比較難實現(xiàn)，因此出現(xiàn)了許多不完全滿足協(xié)調(diào)性要求的“非協(xié)調(diào)單元”或“部分協(xié)調(diào)單元”，有時它們的精度也很好。 插值多項式選擇條件 n插值多項式應該盡可能滿足其收斂性條件（收斂性）n由插值多項式所確定的場函數(shù)變化應該與局部坐標系的選擇無關（各向同性） n假設的插值多項式系數(shù)的數(shù)量應該等于單元的節(jié)點數(shù)（解的唯一性） 選擇條件插值多項式選擇條件n深入分析由收斂性條件可知，插值多項式中必須含有常數(shù)項（剛體位移項），高階項的次數(shù)必須依次增加，不允許有跳躍332210)(xxxxu3310)(xxxu插值多

19、項式選擇條件由選擇條件可知，插值多項式函數(shù)在所有自由度方向上要滿足各向同性性，這樣就不會隨局部坐標系變化而改變了 n深入分析xyyxy, xu3210)(23210)(xyxy, xu插值多項式選擇條件n深入分析選擇條件是為了能由單元節(jié)點值唯一確定插值多項式 4節(jié)點四邊形的插值多項式應該是 xyyxy, xu3210)(插值多項式系數(shù)i (i = 0,1,2,3) 也是4個 單元模型構造整體坐標系法n基本思想 針對彈性體有限元網(wǎng)格建立一個統(tǒng)一的坐標系，每個單元的插值多項式都在這個坐標系上建立 y x o 彈性體 單元模型構造整體坐標系法n2節(jié)點線單元12 oxu1u2x1x2ux1. 假設插值

20、多項式xaaxu10)(2. 利用節(jié)點值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua單元模型構造整體坐標系法3. 代入a0 和 a1，得插值多項式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同類項，設 l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu單元模型構造整體坐標系法lxxN21Nu21212211)(uuNNuNuNxulxxN12N1 和 N2 稱為單元的形函數(shù)；N 稱為單元的形函數(shù)矩陣；u 稱為單元節(jié)點位移向量。 n2節(jié)點線單元的形函數(shù)單元

21、模型構造整體坐標系法n二維3節(jié)點三角形單元 (x1 , y1) (u1 , v1) 1 y x o v 3 2 (x2 , y2) (x3 , y3) u (u3 , v3) (u2 , v2) (u , v) (x , y) 建立整體坐標系oxy 單元模型構造整體坐標系法1. 假設插值多項式y(tǒng)xyxu210),(2. 首先，利用節(jié)點值求 0 、 1 和 2 n二維3節(jié)點三角形單元 323103222102121101yxuyxuyxuyxyxv210),(單元模型構造整體坐標系法)(21)(21)(21332211233221113322110ucucucAubububAuauauaA332

22、21111121yxyxyxA A為單元面積單元模型構造整體坐標系法12332123132ax yx ybyycxx23113231213ax yx ybyycxx31221312321ax yx ybyycxx3. 將 0 、 1 和 2 代入插值多項式，按u1、u2、u3合并同類項332211),(uNuNuNyxu單元模型構造整體坐標系法)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN4. 同理可得332211),(vNvNvNyxv單元模型構造整體坐標系法5. 單元插值多項式為332211332211),(),(vNvNvNyxvuNuNuNy

23、xu)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN單元模型構造整體坐標系法332211321321000000vuvuvuNNNNNNvu6. 單元插值多項式寫成矩陣形式（常用）單元模型構造整體坐標系法321321321321000000vvvuuuNNNNNNvu7. 單元插值多項式的另一種矩陣形式（不常用）單元模型構造整體坐標系法n4節(jié)點四面體單元 2 1 (u1 , v1 , w1) x y z 4 3 (x1 , y1 , z1) (x3 , y3 , z3) (x4 , y4 , z4) (u3 , v3 , w3) (u2 , v2 ,

24、w2) (x2 , y2 , z2) (u4 , v4 , w4) (x , y , z) (u , v , w) 單元模型構造整體坐標系法zyxzyxwzyxzyxvzyxzyxu421042104210),(),(),(1. 假設插值多項式443322114433221144332211),(),(),(wNwNwNwNzyxwvNvNvNvNzyxvuNuNuNuNzyxu2. 插值多項式為單元模型構造整體坐標系法)(61zdycxbaVNiiiii（i=1,2,3,4） 4443332221zyxzyxzyxa 4433221111zyzyzyb4433221111zxzxzxc111

25、4433221yxyxyxd 444333222111111161zyxzyxzyxzyxV 循環(huán)輪換腳標1、2、3、4，相應可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4 單元模型構造整體坐標系法444333222111432143214321000000000000000000000000wvuwvuwvuwvuNNNNNNNNNNNNwvu3. 單元插值多項式寫成矩陣形式（常用）單元模型構造整體坐標系法4. 單元插值多項式另一種矩陣形式（不常用）432143214321432143214321000000000000

26、000000000000wwwwvvvvuuuuNNNNNNNNNNNNwvu單元模型構造整體坐標系法n從理論上講，整體坐標系法可以求任意單元的形函數(shù)，但計算過程太復雜n只能求一維2節(jié)點線單元、二維3節(jié)點三角形單元和三維4節(jié)點四面體單元3種簡單單元的形函數(shù)n復雜的或二次以上的單元必須采用局部坐標系法求n位移場 u 是形函數(shù) Ni 的線性組合，因此形函數(shù)Ni同樣具有插值多項式的特性單元模型構造整體坐標系法n單元形函數(shù)的特性正規(guī)性：單元形函數(shù)之和等于1。 正交性：形函數(shù)在本節(jié)點的值等于1，在其它節(jié)點的值等于0。 單元剛度矩陣2節(jié)點線單元n一維2節(jié)點線單元n單元插值關系 eNuu n單元幾何關系 n

27、單元本構關系 Lu DeN=N1 N2 De=E21uueuxLlxxNlxxN/ )(/ )(1221單元剛度矩陣2節(jié)點線單元n單元剛度矩陣vedvBDBkTvedvBDBkT x yzedxdydzBDBTxedxABDBTA為單元截面積；l為單元長度n矩陣BLNB LNB /21xNxN/1/1llBDBeAlT1111lAE/1/1/1/1llEllAl單元剛度矩陣三角形單元n二維3角形單元n單元插值關系 eNuu T332211vuvuvueu321321000000NNNNNNN)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN單元剛度矩陣三

28、角形單元n單元幾何關系 Lu Luvuxyyxxyyyxx00單元剛度矩陣三角形單元n單元本構關系 DexyyyxxxyyyxxE2100010112平面應力問題單元剛度矩陣三角形單元n矩陣BLNB LNB xNyNxNyNxNyNyNyNyNxNxNxN33221132132100000033221132132100000021bcbcbccccbbbA單元剛度矩陣三角形單元n單元剛度矩陣vedvBDBkTvedvBDBkT x yzedxdydzBDBTx yedxdyhBDBTh為單元厚度BDBehATk為對稱的6*6常數(shù)矩陣A為單元面積單元模型等參單元n等參單元 單元內(nèi)任意一點的位移u

29、與單元節(jié)點位移ue之間的關系為 eNuu 一般單元坐標的插值關系也采用與位移插值關系相同的變換關系即單元內(nèi)任意一點的坐標x與單元節(jié)點坐標xe之間的關系為 eNxx 單元模型等參單元n等參單元凡是幾何形狀和位移場采用同階同參數(shù)插值關系來描述的單元，稱為等參單元 前面介紹的所有單元都屬于等參單元 在描述單元的幾何形狀和位移場時，并不一定非采用同階插值關系 單元模型等參單元332211321321000000vuvuvuNNNNNNvu332211321321000000yxyxyxNNNNNNyxn等參單元3節(jié)點三角形等參單元 單元模型等參單元n超參單元如果幾何形狀插值函數(shù)的階數(shù)高于位移場插值函數(shù)

30、的階數(shù)，稱為超參單元 n次參單元如果幾何形狀插值函數(shù)的階數(shù)低于位移場插值函數(shù)的階數(shù)，稱為次參單元 單元平衡方程組裝過程 n為什么要組裝 ? 2 F 1 3 消除內(nèi)力n組裝的原則是什么 ? 單元自由度與結構自由度對應單元平衡方程組裝過程 2 F 1 3 U3U4U2U1U5U6654321UUUUUUU結構自由度向量U單元平衡方程組裝過程 3 U611 U2U12u1u2U5u3u443214321UUUUuuuueu65214321UUUUuuuueu2 1 U3U4U2U11u1u2u3u42單元平衡方程組裝過程1413121143211441431421411341331321311241

31、2312212111411311211143214321ffffuuuukkkkkkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u42單元平衡方程組裝過程 3 U611 U2U12u1u2U5u3u424232221432124424324224123423323223122422322222121421321221165216521ffffuuuukkkkkkkkkkkkkkkk單元平衡方程組裝過程n單元矩陣擴階11111112131411112122232411113132333411114142434412345610020030040050000006000000kkkk

32、kkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u42單元平衡方程組裝過程 3 U611 U2U12u1u2U5u3u4222211121314222221222324222231323334222241424344 1 2 3 4 5 610020030000004000000500600kkkkkkkkkkkkkkkkn單元矩陣擴階單元平衡方程組裝過程n組裝總剛12121122111112121314131412121122212122222324232411113132333411114142434422223132333422224142434412 3456123004

33、00500600kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk單元平衡方程組裝過程n總體平衡方程1212112211111121213141314121211222212122222324232411113313233341111441424344222231323334522224142434461 2 345612300400500600UkkkkkkkkUkkkkkkkkUkkkkUkkkkkkkkUkkkkU1211122213142324ffffffff單元平衡方程組裝過程00000000000000000000006543216543211413121165432

34、1144143142141134133132131124123122121114113112111ffffUUUUUUkkkkkkkkkkkkkkkk組裝單元單元平衡方程組裝過程242314132212211165432124424324224123423323223114414314214113413313213122422312412322212222112121421311411321211221111100000000654321654321ffffffffUUUUUUkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk再組裝單元FKU 總體剛度方程 K 稱為總體剛度矩陣

35、U 稱為位移向量 F 稱為載荷向量 總體剛度矩陣K的特性 n對稱性 n奇異性 n稀疏性 n非零元素帶狀分布 約束處理過程 n為什么要約束處理 ?總體平衡方程組是奇異的消除無限制的剛體運動 使總體平衡方程組存在唯一一組解約束處理過程邊界條件n 邊界條件分類 力(載荷)邊界條件位移邊界條件 集中載荷力 表面分布力 自重力熱交換引起的溫度載荷 固定位移約束 強制位移約束 關聯(lián)位移約束 002.U 005.U 011.U 514.U CkUU78約束處理過程模型簡化yxUU約束處理過程模型簡化yxUxyU約束處理過程約束方程123456789101112yxU0 . 0321UUU0 . 012963

36、VVVVUUUU121110約束處理過程約束處理方法n位移約束處理方法 賦0賦1法 乘大數(shù)法 約束處理過程賦0賦1法n位移約束處理例子FKU 654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKCkUU34關聯(lián)約束方程 約束處理過程賦0賦1法n賦0賦1法6543216543216665646362615655545352514645444342413635343332312625242

37、32221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK646545444343242141654321666564636261565554535251464544434241363534333231262524232221161514131211000000CKFCKFCKFCKFCKFCKFUUUUUUKKkKKKKKKkKKKKKKkKKKKKKkKKKKKKkKKKKKKkKKKK約束處理過程賦0賦1法n有6個方程，5個未知數(shù)，如果約束方程可以消除有限元平衡方程組的奇異性，則取任意5個方程聯(lián)立求解，都會得到方程組

38、的唯一一組解。 n系數(shù)矩陣由原來的對稱的變成了非對稱的，這對于大規(guī)模有限元方程組求解是十分不利的，采用相同的求解方法，在求解時間和矩陣存貯容量方面都增加了一倍。 約束處理過程賦0賦1法n如何解決 ?可以發(fā)現(xiàn), 約束處理過程實際上是在系數(shù)矩陣上做了一次列初等變換。 為了保持平衡方程組系數(shù)矩陣的對稱性，對變換后的系數(shù)矩陣再做一次相同的行初等變換。 約束處理過程賦0賦1法6465454443432421416543216665646362615655545352514636453544433433423241312625242322211615141312110)(000010000)(00CKFC

39、KFCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkKKkKKkKKkkKKkKKkKKKKkKKKKKKkKKKK具體做法：第4行乘以系數(shù)k加到第3行，并去掉第4行。為了保持系數(shù)矩陣的階數(shù)，將第4行的所有元素賦0，在其對角線元素位置賦1。即所謂賦0賦1法。 646545444444343242141654321666564636261565554535251464544434241463645354443343342324131262524232221161514131211)(0000)(00CKFCKFCKFCKFkCKFCKFCKFUUUUUUKKkKKKKKKkK

40、KKKKKkKKKKkKKkKKkKKkkKKkKKkKKKKkKKKKKKkKKKK約束處理過程賦0賦1法n經(jīng)過初等變換，方程組的系數(shù)矩陣仍然保持對稱性 n初等變換不會改變方程組的解 n約束后的方程組可以求得5個未知數(shù) n通過關聯(lián)約束方程回代求解U4約束處理過程賦0賦1法n進一步改進方法 將關聯(lián)約束方程直接引進約束后的方程組中，使U4也直接參與方程組求解。 用關聯(lián)約束取代約束后方程組的賦0賦1行（第4行），然后再做對稱化處理。即將取代后的一行方程（第4行）乘以k加到第3行，則系數(shù)矩陣仍然保持對稱性。約束處理過程賦0賦1法646545444343242141654321666564636261

41、5655545352514636453524443343342324131262524232221161514131211)(0000100)(00CKFCKFCkCCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkkKKkKKkkkKKkkKKkKKkKKKKkKKKKKKkKKKK約束處理過程賦0賦1法n基本原理 利用初等變換對求解方程組進行相同的行列變換，既保證方程組解不會改變，又可以保持方程組系數(shù)矩陣的對稱性。 在進行初等變換時，只要保證對方程組系數(shù)矩陣做相同的行列變換，就可以保持方程組系數(shù)矩陣的對稱性。 約束處理過程賦0賦1法n固定和強制位移約束條件處理 通過關聯(lián)位

42、移約束方法簡化。如果k = 0，則退化成強制位移邊界條件即 U4=C 如果k =C= 0，則退化成固定位移邊界條件即 U4=0 6465454443432421416543216665646362615655545352514636453524443343342324131262524232221161514131211)(0000100)(00CKFCKFCkCCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkkKKkKKkkkKKkkKKkKKkKKKKkKKKKKKkKKKK約束處理過程賦0賦1法n強制位移約束條件處理 U4=C646545343242141654321

43、6665636261565553525136353332312625232221161513121100001000000CKFCKFCCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK約束處理過程賦0賦1法n固定位移約束條件處理 U4=06543216665636261565553525136353332312625232221161513121100001000000UUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK653210FFFFF646545343242141CKFCKFCCKFCKFCKF約束處理過程賦0賦1法n固定位移和強制位移約束處理

44、后的系數(shù)矩陣是相同的，只是簡單地將方程組系數(shù)矩陣中要約束自由度的行列分別賦0，對角線元素賦1，這也是賦0賦1法的由來。約束處理過程賦0賦1法n在方程組載荷右端項的處理方法上兩者是不同，處理固定位移邊界條件時，只要將對應自由度的載荷賦0即可。處理強制位移邊界條件時，要在方程組系數(shù)矩陣未賦0賦1法前，先將對應自由度的列乘以系數(shù)C減到載荷右端項，再將對應自由度的載荷位置賦C。約束處理過程乘大數(shù)法n乘大數(shù)法n基本原理 利用矩陣的初等變換不改變方程組解的思想。 約束處理過程乘大數(shù)法65432165432166656463626156555453525146454443424136353433323126

45、2524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKCkUU34關聯(lián)約束方程 約束處理過程乘大數(shù)法CUkU43關聯(lián)約束方程 ACAUAkU43A是一個大數(shù)，是系數(shù)矩陣中對角線元素K44的1010倍量級以上 為什么要乘以大數(shù)A ？放大位移約束方程的優(yōu)勢約束處理過程乘大數(shù)法關聯(lián)約束方程 ACAUAkU43將關聯(lián)約束加到有限元平衡方程對應自由度行，第3行或第4行，這里取第4行6543216543216665646362615655545352514645444342413635343332312625242322

46、21161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKKKKKKKKKKKKKKKKK約束處理過程乘大數(shù)法n如果關聯(lián)約束方程可以消除有限元平衡方程的奇異性，約束后的方程組就存在唯一的一組解。n約束后的方程組的系數(shù)矩陣是非對稱的。n利用初等變換方法將系數(shù)矩陣變換成對稱的 約束處理過程乘大數(shù)法關聯(lián)約束方程 ACAUAkU43再乘以系數(shù)-k AkCAkUUAk432加到約束后方程組的第3行 約束處理過程乘大數(shù)法654321654321666564636261565554535251464544434241363534333231262524232221

47、161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKKKKKKKKKKKKKKKKK6543216543216665646362615655545352514645444342413635342333231262524232221161514131211FFACFAkCFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKAkKAkKKKKKKKKKKKKKKK約束處理過程乘大數(shù)法n強制位移邊界條件 )0(4kCU約束后的方程組簡化為 6543216543216665646362615655545352514645444342413635

48、34333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK6543216543216665646362615655545352514645444342413635342333231262524232221161514131211FFACFAkCFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKAkKAkKKKKKKKKKKKKKKK約束處理過程乘大數(shù)法n固定位移邊界條件 k = 0，C = 0 約束后的方程組簡化為 65432165432166656463626156555

49、4535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK約束處理過程乘大數(shù)法n固定位移和強制位移邊界條件的乘大數(shù)約束處理相對比較簡單，而且它們的系數(shù)矩陣約束后是相同的，只是簡單地

50、將方程組系數(shù)矩陣中要約束自由度的對角線元素加上一個相對大數(shù)A即可 n乘大數(shù)法的叫法并不十分準確，應該叫加大數(shù)法更貼切 n乘大數(shù)和加大數(shù)的效果是一樣的約束處理過程兩種方法比較n賦0賦1法在約束處理過程中是嚴格精確的，而乘大數(shù)法是一種近似約束處理方法，它的精度取決于所乘大數(shù)A值 兩種方法都可以消除有限元平衡方程的奇異性，得到符合實際邊界條件的唯一一組解。但兩種方法還是有很大的區(qū)別 約束處理過程兩種方法比較n采用乘大數(shù)法約束處理后的有限元平衡方程在求解時可能造成解的失真，大數(shù)A值越大可能解的偏差會越大，而賦0賦1法就不會出現(xiàn)類似的問題，它在約束過程和求解過程都是精確的n乘大數(shù)法相對于賦0賦1法在約束

51、處理過程上簡單一些 約束處理過程兩種方法比較n賦0賦1法實際上是將關聯(lián)位移約束方程代入到有限元平衡方程中的，是代入法。而乘大數(shù)是將占絕對優(yōu)勢的關聯(lián)位移約束方程合并到有限元平衡方程中的，是罰方法，計算誤差來自于合并過程，計算精度取決于關聯(lián)位移約束方程的優(yōu)勢大小n商業(yè)軟件中，位移邊界條件的約束處理都采用賦0賦1法，乘大數(shù)很少被采用主要原因是它是一種近似方法，而且大數(shù)的大小也不好確定，有時還會造成求解失敗 約束處理過程彈簧單元假設柔性彈簧 kOXYU4f f = kU4k約束處理過程彈簧單元彈簧約束方程 f = kU465432165432166656463626156555453525146454

52、4434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFfFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK6544321654321666564636261565554535251464544434241363534333231262524232221161

53、514131211FFkUFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKkKKKKKKKKKKKKKKKKKKKKKK方程組求解過程特點n方程組求解是有限元計算過程中很重要的一部分，在有限元法的發(fā)展過程中，有限元方程的求解效率一直是其應用的最大瓶頸之一 n有限元方程組的特點： 有限元方程組的系數(shù)矩陣具有對稱、稀疏、帶狀分布以及正定、

54、主元占優(yōu)。有效地利用這些特點，以減少系數(shù)矩陣的存貯量，提高方程組求解效率 方程組求解過程分類比較n線性方程組的解法主要分兩大類: 直接解法：以高斯消去法為基礎，以等帶寬或變帶寬方式存貯系數(shù)矩陣內(nèi)元素，對于求解規(guī)模比較大的問題，要存貯的元素非常巨大。 迭代解法：只需要存貯系數(shù)矩陣中非零元素，存貯量很小，一般是變帶寬存貯量的20%或更少，有些算法的求解效率也非常高，適合求解大規(guī)模線性方程組。但是這種解法對接近病態(tài)的方程組很難保證收斂性。 方程組求解過程帶寬定義n有限元方程組系數(shù)矩陣是稀疏的、非零元素呈帶狀分布，帶寬就是它的寬度，帶寬的大小是由系統(tǒng)有限元網(wǎng)格的節(jié)點號排序決定的，具體求法是 帶寬=(單

55、元最大節(jié)點號之差+1)*節(jié)點自由度數(shù) n帶寬是網(wǎng)格節(jié)點標注方法直接決定的，不同標注方法帶寬可能相差很大 方程組求解過程帶寬n帶寬是網(wǎng)格節(jié)點標注方法直接決定的，不同標注方法帶寬可能相差很大 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22

56、 23 24 28 26 27 25 方程組求解過程帶寬n所示四邊形網(wǎng)格的三種節(jié)點號標注方法，每個節(jié)點是2個自由度n結構的帶寬分別是12，18，56，相差很大，其中12和56之間相差近5倍，這就意味著系數(shù)矩陣的存貯量也是相差5倍，因此，對于大規(guī)模復雜系統(tǒng)的節(jié)點號優(yōu)化是十分必要的 方程組求解過程系數(shù)矩陣存貯 n系數(shù)矩陣存貯 如果節(jié)點號排序優(yōu)化的比較好，系數(shù)矩陣的存貯量就會減少很多。根據(jù)系數(shù)矩陣的對稱性，一般都是按半帶寬存貯。n系數(shù)矩陣存貯的方法 二維等帶寬存貯 一維變帶寬存貯 方程組求解過程二維等帶寬存貯 n二維等帶寬存貯 D D 0 0 K K n 方程組求解過程二維等帶寬存貯n二維等帶寬存貯

57、消除了最大帶寬以外的全部零元素，節(jié)省了系數(shù)矩陣元素的存貯量。但是由于取最大帶寬為存貯范圍，因此不能排除在帶寬內(nèi)的大量零元素。當系數(shù)矩陣的各行帶寬變化不大時，適合采用二維等帶寬存貯，方程組求解過程中系數(shù)矩陣元素的尋址也比較方便，求解效率較高。n當出現(xiàn)局部帶寬特別大的情況時，采用二維等帶寬存貯時，將由于局部帶寬過大而使整體系數(shù)矩陣的存貯大大增加。 方程組求解過程一維變帶寬存貯 n一維變帶寬存貯 一維變帶寬存貯方法就是把變化的帶寬內(nèi)的元素按一定的順序存貯在一個一維數(shù)組中。由于它不按最大帶寬存貯，因此比二維等帶寬存貯更節(jié)省內(nèi)存。按照解法可分為按行一維變帶寬存貯和按列一維變帶寬存貯。 n按行一維變帶寬存

58、貯 對稱 20181296431M方程組求解過程一維變帶寬存貯 輔助的尋址數(shù)組M n一維變帶寬存貯是最節(jié)省內(nèi)存的一種方法，但是由于要借助于尋址數(shù)組尋找系數(shù)矩陣元素的位置，相對二維等帶寬存貯方法來說要復雜一些，而且在程序?qū)崿F(xiàn)時也要復雜得多，方程組求解過程中也要消耗一些數(shù)組尋址時間。因此，在選用存貯方法時要權衡二者的利弊，統(tǒng)盤考慮。一般當帶寬變化不大，計算機內(nèi)存允許時，采用二維等帶寬存貯方法是比較合適的。 方程組求解過程一維變帶寬存貯 方程組求解過程求解方法n方程組求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-賽德爾（Gauss-Seidel）迭代法 應變、應力回代過程 n

59、單元應變和應力回代求解 通過求解有限元平衡方程得到有限元節(jié)點位移后，就可以進行系統(tǒng)的剛度校核。如果所分析問題要進行強度校核，就要回代求解單元的應變和應力。n由插值關系和幾何關系可得單元應變，再通過本構關系得到單元應力eeBuLNuDe板料成形有限元法分類n彈塑性有限元法n剛塑性有限元法彈塑性剛塑性板料成形有限元法分類n靜力隱式有限元法n靜力顯式有限元法n動力顯式有限元法n逆算法(一步成形)板料成形有限元法單元模型n單元模型n 薄膜單元n 薄殼單元n 中厚殼單元n 等效彎曲單元板料成形有限元法單元平衡方程n單元平衡方程n單元插值關系(局部坐標系) eNuu n單元幾何關系(局部坐標系) n單元本

60、構關系(局部坐標系) Lu De0TTTvavPdvdadvGuPun最小勢能原理(局部坐標系) n坐標變換關系 eeTUu 板料成形有限元法單元平衡方程0TTTvavPdvdadvGuPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu代入局部坐標系的單元3個關系式, 得代入坐標轉(zhuǎn)換關系式, 得0)()()(TTTTTTTTTvevaeeeedvdadvGNTUPNTUTUBDBTU0TTTTTTvvaeedvdadvGNTPNTTUBDBT0)()()(TTTTTTvvaeedvdadvGNTPNTTUBDBTeeefTTUkTTTeeeFUK板料成形有限元法單