有限元法的前后處理

《有限元法的前后處理》由會員分享，可在線閱讀，更多相關(guān)《有限元法的前后處理（46頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、8 有限元法的前后處理8-1 引言有限元法分成以下幾步：（1）計算模型的幾何剖分、數(shù)據(jù)生成和數(shù)據(jù)準(zhǔn)備；（2）計算分析；（3）計算結(jié)果的分析、整理和圖形顯示。有限元發(fā)展初期，（1）、（3）都是由人工來做。有統(tǒng)計資料表明，（1）的工作量為45%，（3）的工作量為50%，而（2）只占到5%。 以（1）為例，首先按比例畫到坐標(biāo)紙上，然后按一定的順序編號，再整理出如下信息，以供有限元分析之用：*節(jié)點信息節(jié)點編號和節(jié)點坐標(biāo)；*單元信息單元編號和單元中節(jié)點號的排列順序；*材料信息計算模型的材料性質(zhì)（彈性模量，泊松比，比熱，導(dǎo)熱系數(shù)等）；*載荷信息計算模型所受的負(fù)載信息（集中力，體積力，表面力，溫度，壓力等）

2、；*約束信息初始條件和邊界條件。 最后還要將這些數(shù)據(jù)一一輸入計算機。以上工作乏味而且容易出錯。有限元分析的前處理就是使計算機部分或全部完成計算模型的幾何剖分、數(shù)據(jù)生成和數(shù)據(jù)輸入，有限元的后處理則是將有限元計算結(jié)果由計算機整理成易于閱讀或分析的數(shù)值或圖形形式。8-2 有限元分析的前處理技術(shù)有限元法的前處理主要有以下內(nèi)容：（1）計算模型的幾何表示；（2）模型網(wǎng)格的自動分劃（或剖分）；（3）剛度矩陣的帶寬優(yōu)化；（4）模型網(wǎng)格圖的計算機繪制。 一 計算模型的幾何表示1 對于設(shè)計階段的零部件進行有限元分析，如果采用計算機輔助設(shè)計，則模型的幾何表示可采用計算機造型系統(tǒng)中的幾何表示。常用的有表面模型和立體模

3、型。2 對于已有的零部件做有限元分析，既可以將其輸入到計算機中，采用第一種辦法表示其幾何形狀，也可以根據(jù)零件的幾何形狀來決定其表示法。常用的有整體表示法和分塊表示法。（1）整體表示法：用點表示線，用線表示面，用面表示體。（2）分塊表示法：把整體看成是由簡單個體的組合，而簡單個體則可用多邊形和多面體表示。 二 有限元網(wǎng)格的自動剖分有限元網(wǎng)格的自動剖分與計算模型的幾何表示方法有密切的關(guān)系。整體表示的幾何模型，適宜于采用整體剖分要用到較多的數(shù)學(xué)知識。下面討論分塊表示幾何模型時常用的一種分塊剖分法。以二維問題為例。圖示的二維區(qū)域可以看成是A,B,C三個部分組成。 A,B,C都是簡單的四邊形，四邊形頂點

4、的坐標(biāo)可以表示其形狀和位置。 網(wǎng)格自動剖分的方法分成四個步驟：1）分塊映射；2）網(wǎng)格剖分；3）順序編碼；4）總體合成。1 分塊映射首先把計算區(qū)域手工粗分成若干個四邊形區(qū)域的組合，每個四邊形稱為一個大單元。 B CA 為了既能表示直邊大單元，又能表示曲邊大單元，采用8個節(jié)點來描述一個大單元。右圖中，（a）是真實圖形在總體坐標(biāo)系下的樣子，經(jīng)過等參數(shù)單元的變換后，得到（b）圖中邊長為2的正方形。同樣的方法可以得到B和C的映射結(jié)果。 2 3 467 11 2 3 4 5567 88yO OA )(b)(a x 2 網(wǎng)格剖分網(wǎng)格剖分是在映射后的正方形區(qū)域中進行的。剖分時可以有兩種方法：等邊剖分和不等邊剖

5、分。（1）等邊剖分將=-1邊和=1邊剖分成相等的份數(shù)。同樣，將=-1和=1邊也剖分成相等的份數(shù)。則最終剖份單元數(shù)為 。 NN O （2）不等邊剖分某一方向兩個邊界上的剖分?jǐn)?shù)不一致。比如 =-1分割數(shù)為3， =1一側(cè)的分割數(shù)為7。需要計算沿另一個方向（等值、 方向）的分割數(shù)，使相鄰的兩條等值線的分割數(shù)差1。 O3 4 6 75CBA 對于相鄰的兩條等值線分割數(shù)不相等的情況，應(yīng)將分割數(shù)多的等值線上的最后兩點與相鄰的、分割數(shù)少1的等值線上的最后一點組成一個三角形單元。已經(jīng)形成的四邊形單元，只要連接四邊形的兩個對角點，即可得到兩個三角形單元?；蛘弑３炙倪呅螁卧@樣，求解區(qū)域既有四邊形單元，也有三角形

6、單元，相應(yīng)的單元剛度矩陣也有兩種類型，總體合成時，按照總體節(jié)點編號累加單元中的元素值即可。完成局部坐標(biāo)系o中的大單元的分割后，對應(yīng)的總體坐標(biāo)系xoy中的大四邊形也同時被分割了。從o到xoy的坐標(biāo)變換關(guān)系是 式中 為映射函數(shù)，也就是8節(jié)點四邊形等參數(shù)單元位移插值函數(shù)的形函數(shù)。 是總體坐標(biāo)系中大四邊形各邊節(jié)點的坐標(biāo)。 ii ii ii yNy xNx 8181 , 8,2,1,iNi ,84 1121 2,6 1121 7,5,3,1 11141, 22 iiiN ii iiiii 8,2,1,iyx ii 按照坐標(biāo)變換關(guān)系式，將已經(jīng)分割后的局部坐標(biāo)系中的正方形中的網(wǎng)格節(jié)點坐標(biāo)(,)代入，即可得到

7、對應(yīng)的總體坐標(biāo)系下的節(jié)點坐標(biāo)(x,y)。例如， 方向分割數(shù)N=4， 方向的分割數(shù)N=3，求圖中K點的坐標(biāo)(xK,yK)。 KO 如果以O(shè)為起點，則K點的局部坐標(biāo)為根據(jù)已知的 ， 將(K,K)代入坐標(biāo)變換式，即可求出對應(yīng)的(x K,yK)。 方向的分割線序數(shù)。方向和點算起從，；起算點的局部坐標(biāo) OKK KN KNOO OK OK 1,1, 3113212 2 114212 8,2,1,iN i 8,2,1,iyx ii 3 順序編碼編碼就是將所有節(jié)點連續(xù)無缺地順序編號，并對每個單元編號和形成單元信息。編碼可以分成兩步來做。首先對每個大四邊形編碼。對基本的大四邊形單元編碼時，可以按行（列）進行。以

8、按行編碼為例，每計算一個節(jié)點，節(jié)點號加1。單元編號也可類似進行（比如按列）。 )2( )3(1k2k )1( )4( )5(5k )1(NS )2(NS )3(NS )4(NS 為了形成與某一個單元相關(guān)的節(jié)點號，在節(jié)點編號時，對每一個新行的第一個節(jié)點專門形成一個數(shù)組NS(i)。其中的i代表節(jié)點編碼的行號， NS(i)的值代表該行起點節(jié)點號。如果單元編碼按列進行，則第二個單元的節(jié)點是NS(3)， NS(3)+1， NS(2)+1， NS(2)。如果單元是8節(jié)點等參元，還要計入單元每邊中間一點的節(jié)點號。對每個大四邊形順序編碼，這樣，可能會有一些節(jié)點在大四邊形組合后，同一節(jié)點擁有兩個節(jié)點號。這個問題

9、可以在下一步總體合成中予以解決。 4 總體合成將經(jīng)過剖分、編碼后的大四邊形單元合在一起，調(diào)整具有多個節(jié)點號的節(jié)點編碼，重排節(jié)點號，稱為總體合成?？傮w合成的主要任務(wù)是使整體節(jié)點號的編號連續(xù)無缺，并修改分開編碼形成的單元信息。對于同一節(jié)點具有多個節(jié)點編號多重編號的節(jié)點，可以采用比較節(jié)點坐標(biāo)的辦法將其消除。比如，如果同時成立，說明i和j重合或很接近。 yji xji yy xx 21 21 式中(xi1,yi1)為第一大單元中的節(jié)點i的坐標(biāo)；(xj1,yj1)為第二大單元中的節(jié)點j的坐標(biāo)； x,y是預(yù)先給定的距離很小的數(shù)值。找出重復(fù)編號的節(jié)點之后，修改第二大單元中的節(jié)點信息和單元信息：（1）保留節(jié)點

10、號i，取消節(jié)點編號j；（2）將第二大單元中比i大的節(jié)點編號全部減1；（3）修改第二大單元中該節(jié)點所在的單元信息和比該節(jié)點編號要大的那些單元的信息。重復(fù)以上過程，直至第二大單元中的所有節(jié)點都循環(huán)完畢。 對于三維問題，也可以實現(xiàn)有限元網(wǎng)格的自動分割，只是分塊區(qū)域為六面體或五面體，映射后的區(qū)域為正六面體和直三棱柱。 1 2 3 4 56 1016 1511 17 18 19 2021 2224 2627 2930 3233 34 35 1823262932 33 2528大單元連續(xù)編號重復(fù)編號的剔除222123 3134 三 剛度矩陣的帶寬優(yōu)化為了節(jié)省計算機的存儲空間，需要優(yōu)化總剛矩陣的帶寬。節(jié)點編

11、號對應(yīng)總剛矩陣中的非零元素的位置。節(jié)點編號變動，總剛矩中非零元素的位置也發(fā)生變動，反之亦然。帶寬優(yōu)化的原理：通過調(diào)整總剛矩陣中非零元素的位置，對應(yīng)地修改單元信息，從而減少剛陣帶寬。 為此引入鄰接矩陣B，其階數(shù)與剛度矩陣相同，其中的元素非0即1式中，aij是剛度矩陣中的元素。B也具有帶狀的樣子。帶寬優(yōu)化有許多實用算法，有些要用到圖論或較深的數(shù)學(xué)知識。下面介紹兩種易于理解的帶寬優(yōu)化方法，它們都是采用變換鄰接矩陣中的行或列的辦法來減少剛陣帶寬的。 0 0 0 1 ijijij aab 1 羅森算法羅森（Roson）算法的基本思想：找出引起B(yǎng)矩陣最大帶寬的一對頂點，取其中一個與其他頂點互換，看能否減小

12、帶寬。若可以，則進行交換。交換后，再繼續(xù)尋找引起當(dāng)時B矩陣最大帶寬的一對頂點，并重復(fù)上述過程。為此，首先將B中的“頂點”序號紀(jì)錄在一個數(shù)組中，稱為“頂點表”?！绊旤c”是B中每一行離主對角線元素最遠的那一點。然后在頂點表中：（1）確定B矩陣的帶寬及兩個達到最大帶寬的第一對頂點。在這對頂點中，若較高編號頂點能與一低編號頂點互換以減小帶寬，則轉(zhuǎn)向（6）； （2）若較低編號頂點能與一高編號頂點互換以減小帶寬，則轉(zhuǎn)向（6）；（3）若較高編號頂點能與一低編號頂點互換二保持帶寬不變，則轉(zhuǎn)向（6）；（4）若較低編號頂點能與一高編號頂點互換而保持帶寬不變，則轉(zhuǎn)向（6）；（5）如果在（3）、（4）步中已連續(xù)執(zhí)行過

13、規(guī)定的經(jīng)驗次數(shù)，或已交換過的頂點又重新選出來交換，則算法終止；否則，轉(zhuǎn)向（6）；（6）執(zhí)行所指出的頂點交換，即分別交換頂點表的對應(yīng)分量以及B矩陣的行和列，并轉(zhuǎn)向（1）。可見，整個交換過程都在頂點表中進行。 2 阿基茨厄特庫算法阿基茨(Akyuz)厄特庫(Utku)算法也在1968年提出一種減小平均帶寬的方法，簡稱AU算法。平均帶寬的定義：AU算法的基本步驟：（1）從B矩陣中取出相鄰的兩行列進行交換，并計算平均帶寬。如果滿足下列兩個條件之一，則交換有效：行的帶寬。矩陣第剛度矩陣的階數(shù)；式中in n ini i 11 平均帶寬減少；平均帶寬保持不變，但有較多元素的行從矩陣中心向外移。（2）在一個指

14、定的交換循環(huán)內(nèi)，按(1) 及(1) 執(zhí)行行列交換，其交換順序規(guī)定為(1,2), (n,n-1), (2,3), (n-1,n-2),直至中心行。（3）如果在一個循環(huán)內(nèi)沒有發(fā)生交換，或經(jīng)過經(jīng)驗次數(shù) 次的循環(huán)而平均帶寬不減小，交換運算停止；否則，重復(fù)執(zhí)行(1)、(2)。一般來講，實用的帶寬優(yōu)化算法并不能使剛度矩陣的帶寬達到極小化，但并不影響其實用性。大型的有限元分析軟件包一般都具有帶寬極小化的功能。1003 n 四 模型網(wǎng)格的自動生成目的：檢查所生成的幾何數(shù)據(jù)的正確性。方法：根據(jù)單元數(shù)組信息（單元和單元節(jié)點關(guān)系）和節(jié)點數(shù)組信息（節(jié)點號和節(jié)點坐標(biāo)關(guān)系），將相鄰節(jié)點連接起來，一個單元一個單元地去形成有

15、限元剖分網(wǎng)格。注意：所得結(jié)果還要消除隱藏線，否則看不清結(jié)果。這部分內(nèi)容屬于計算機圖形學(xué)，在此不做討論。 8-3 有限元分析的后處理有限元數(shù)據(jù)后處理的原因：1）計算結(jié)果是位移和應(yīng)力；2）數(shù)據(jù)量大，不直觀。有限元數(shù)據(jù)處理工作包括兩個方面；1）數(shù)值處理 將有限元結(jié)果轉(zhuǎn)化成工程中常用的形式，或設(shè)計師熟悉的形式；2 ）圖形處理 用圖形直觀地表示設(shè)計結(jié)果，使設(shè)計結(jié)果一目了然。 一 數(shù)值處理以應(yīng)力分析為例。有限元位移法計算得到節(jié)點的位移，進而計算單元的應(yīng)力：單元內(nèi)任意一點處的應(yīng)力為 。上述結(jié)果不論是數(shù)值積分點處的應(yīng)力，節(jié)點處的應(yīng)力，還是單元中心處的應(yīng)力，也不論其精度高低，都不好使用，原因是：我們需要的不是應(yīng)

16、力分量，而是主應(yīng)力或等效應(yīng)力。說明：單元內(nèi)不同點處的應(yīng)力值的精度是有差別的。一般，數(shù)值積分點處的應(yīng)力值較為準(zhǔn)確。 )(, eBDzyx zxyzxyzyx , 1 主應(yīng)力計算和等效應(yīng)力計算工程中往往用一點的主應(yīng)力 來計算或判斷結(jié)構(gòu)的安全程度。對于二維問題，有限元分析可得任意一點的 ，則主應(yīng)力 與它們的關(guān)系為主應(yīng)力方向 321 , xyyx 和, 21, xyyxyx xyyxyx 22 21 22 22 yx xy 2arctan21 最大剪應(yīng)力和最小剪應(yīng)力為對于三維應(yīng)力問題，主應(yīng)力可以通過求解下述三次方程得到：2 21minmax 2223 22221 322132 0 xyzzxyyzxz

17、xyzxyzyx zxyzxyxzzyyx zyxIII III 其中 上式的解可表示為將 按從大到小的順序排列則 是最大主應(yīng)力， 是最小主應(yīng)力。 23221 32131221 13 12 11 32 2792cos ,332 34cos3 32cos3 3cos3 II IIIIIIR RI RI RI 其中 321 , 321 1 3 從理論上可以證明，三個根都是實根，因此有否則，一定是有限元計算結(jié)果有誤。主應(yīng)力的方向可由以下三式中的任意兩式三式聯(lián)立，得出主應(yīng)力的方向余弦1cos1- ,03 221 II 1 000 222 iii iiziyzizx iyziiyixy izxixyii

18、x nml nml nml nml以及 3,2,1, inml iii 上式中前兩式與第四式聯(lián)立求解的結(jié)果為令 ，即可得到 的方向余弦。 21222 212 212 12 12 1 11 xyiyixixyzxyzx iyzxyzxy xyiyixi ixyzxyzxi iyzxyzxyiA An Am Al 式中3,2,1i 321 , 另外，要核算所分析的零件強度是否足夠，還要根據(jù)一定的強度理論來計算某點的綜合應(yīng)力。各種強度理論的等效應(yīng)力如下： 213232221xd3 31xd3 321xd2 11xd21 第四強度理論第三強度理論第二強度理論第一強度理論 2 應(yīng)力修勻有限元保證位移在單

19、元邊界的連續(xù)性。但是，應(yīng)力在單元邊界并不連續(xù)。怎樣處理計算結(jié)果才能使應(yīng)力連續(xù)和直觀呢？（1）簡單平均或加權(quán)平均以三角形單元為例。三角形為常應(yīng)力單元，但不同三角形單元的應(yīng)力值不同。簡單平均：加權(quán)平均： 有關(guān)單元面積所在單元面積節(jié)點應(yīng)力值節(jié)點應(yīng)力值有關(guān)單元節(jié)點應(yīng)力值節(jié)點應(yīng)力值 （2）最小二乘法在單元內(nèi)采用最小二乘法修勻應(yīng)力；單元節(jié)點的應(yīng)力值，取環(huán)繞該節(jié)點的不同單元的應(yīng)力值的平均（簡單平均或加權(quán)平均）。單元內(nèi)全應(yīng)力修勻設(shè)應(yīng)力向量 ，修勻后的應(yīng)力向量為 ，則使A(e)達到最小。其中 )( 21)( eV Te dVA 用的插值函數(shù)。值函數(shù)和修勻應(yīng)力時采有限元求解時的單元插節(jié)點應(yīng)力值；應(yīng)力值和修勻應(yīng)力后

20、的有限元求解得到的節(jié)點ii ii iiiiNN NN , , 根據(jù)上式是以修勻后的節(jié)點應(yīng)力向量 為未知量的線性代數(shù)方程組。利用解析法或高斯積分，可以求得 ，進而利用單元應(yīng)力插值，可以計算單元內(nèi)任何一點的應(yīng)力值。 )( )()( ,2,1 0 ,2,1 0 )( eV iT eie nidVN niAe 或 i i 單元內(nèi)部分應(yīng)力修勻只修勻單元中我們感興趣的分量。例如，只修勻分量。這時有x )( 2)( 21eV xxxxe dVA )()( ,2,1 0 )( eV ixxxieii xixi xiixxiix nidVNANN NN e 帶入變分后的表達式函數(shù)。修勻應(yīng)力時采用的插值值函數(shù)；有

21、限元求解時的單元插值；修勻應(yīng)力后的節(jié)點應(yīng)力應(yīng)力值；有限元求解得到的節(jié)點將 上式為一線性代數(shù)方程組，方程的數(shù)目為n(e)，未知量為修勻后的單元節(jié)點應(yīng)力分量 。單元內(nèi)應(yīng)力線性外推利用單元內(nèi)的數(shù)值積分點處的應(yīng)力值，外推單元節(jié)點的應(yīng)力值，然后再平均與節(jié)點有關(guān)的單元節(jié)點應(yīng)力。對于二維四節(jié)點等參數(shù)單元，若計算單元剛度矩陣時取了 個積分點，則積分點的坐標(biāo)為xi22 ,31,31: ,31,31: ,31,31: ,31,31: DC BA 帶入公式可以求出未修勻時積分點處的應(yīng)力值設(shè)修勻后的應(yīng)力在單元內(nèi)按雙線性變化，則修勻應(yīng)力為 )(, eBDzyx DCBA , i i ii NNNN NNNN N 432

22、1 43214321 41, O1 24 3A BCD 令積分點處的應(yīng)力值等于修勻后的應(yīng)力值，則 Pi i DCBADCBAP N NDNDNDNDN CNCNCNCN BNBNBNBN ANANANAN 1 43214321 4321 4321 4321 由此可以解出： 將上述求得的節(jié)點應(yīng)力，再根據(jù)不同的單元節(jié)點應(yīng)力值予以平均，就可得到最終的節(jié)點應(yīng)力值。一般情況下，采用局部應(yīng)力磨平處理，可以得到較好的結(jié)果，而計算量是很小的，所以得到廣泛的應(yīng)用。應(yīng)力全域磨平和子域局部應(yīng)力磨平可參見清華大學(xué)王勖成編著的“有限單元法基本原理和數(shù)值方法”。二 圖形處理圖形處理使計算結(jié)果簡潔，直觀，生動。常用以下幾種

23、方法。 1 變形圖分為兩種：1）用網(wǎng)格表示內(nèi)部變形（當(dāng)外部變形不大時）；2 ）用邊界表示外部變形（只關(guān)注外部變形時）。兩種方法，都把原圖和變形后的圖重疊畫出，以作對比。繪圖時注意：要放大，否則看不出來。例如，原節(jié)點坐標(biāo)(x,y)，位移為(u,v)，則繪制變形時該節(jié)點的坐標(biāo)為 變形放大倍數(shù)。式中QP uQyy uPxx, 2 等值線圖等值線可以是等溫度線，熱流線，等應(yīng)力線。方法：追蹤法+插值法。思路：先找到一個等值點，然后根據(jù)離散節(jié)點坐標(biāo)和節(jié)點的場變量值，求出節(jié)點連線間的等值點位置，并和開始找到的等值點相連。就這樣一直尋找下去，直到邊界。 jA Bia b 入口 入口 出口c 出口 插值算法：設(shè)

24、兩個節(jié)點i和j的坐標(biāo)為(xi,yi),(xj,yj)，場變量(Ti,Tj)，待求場變量為Ta，則當(dāng)成立時，在 邊上有一等值點a，其坐標(biāo)為對于二次插值單元，單元邊界上有三個節(jié)點i,j,k，可采用二次插值或分段線性插值求出等值點。 0 jaia TTTTij ijij iaia ijij iaia yyTT TTyy xxTT TTxx 追蹤算法：從A單元a點到b點，再從b點進入B單元，c為B單元的出口。求解關(guān)鍵：找出單元B。為此，只要搜索除A單元外的其他單元，這些單元的特點是：單元包含節(jié)點i,j；或單元節(jié)點坐標(biāo)與i,j一致。使用中的幾個問題：1）尋找起點沒有規(guī)律可循。按單元循環(huán)，每個單元按邊循環(huán)。2）等值線可能不只一條；3）只搜索沒有搜索過的單元。 3 濃淡色彩圖（彩色云圖）對顯示的區(qū)域不斷剖分，直到一個像素的大小，同時用插值算法計算剖分后的小單元（像素）對應(yīng)的場變量值，然后用當(dāng)前場變量值對應(yīng)的色彩，顯示該像素。4 節(jié)點變量變化圖采用曲線擬合或曲線插值技術(shù)，將所有離散時間點上的某節(jié)點場變量值，按一定的函數(shù)形式擬合成一條或幾條連續(xù)曲線，并按比例畫出。