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è)計階段的零部件進(jìn)行有限元分析，如果采用計算機輔助設(shè)計，則模型的幾何表示可采用計算機造型系統(tǒng)

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

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

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

6、按照總體節(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)。,按照坐標(biāo)變換關(guān)系式，將已經(jīng)分割后的局部坐標(biāo)系中的正方形中的網(wǎng)格節(jié)點坐標(biāo)(,)代入，即可得到對應(yīng)的總體坐標(biāo)系下的節(jié)點坐標(biāo)(x,y)。 例如， 方向分割數(shù)N=4， 方向的分割數(shù)N=3，求圖中K點的坐標(biāo)(xK,yK)。,如果以O(shè)為起點，則K點的局部坐標(biāo)為 根據(jù)已知的 ， 將(K,K)代入坐標(biāo)變換式，即可求出對應(yīng)的(xK,yK)。,3 順

7、序編碼 編碼就是將所有節(jié)點連續(xù)無缺地順序編號，并對每個單元編號和形成單元信息。 編碼可以分成兩步來做。首先對每個大四邊形編碼。 對基本的大四邊形單元編碼時，可以按行（列）進(jìn)行。以按行編碼為例，每計算一個節(jié)點，節(jié)點號加1。單元編號也可類似進(jìn)行（比如按列）。,為了形成與某一個單元相關(guān)的節(jié)點號，在節(jié)點編號時，對每一個新行的第一個節(jié)點專門形成一個數(shù)組NS(i)。其中的i代表節(jié)點編碼的行號， NS(i)的值代表該行起點節(jié)點號。 如果單元編碼按列進(jìn)行，則第二個單元的節(jié)點是NS(3)， NS(3)+1， NS(2)+1， NS(2)。 如果單元是8節(jié)點等參元，還要計入單元每邊中間一點的節(jié)點號。 對每個大四邊

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

9、之后，修改第二大單元中的節(jié)點信息和單元信息： （1）保留節(jié)點號i，取消節(jié)點編號j； （2）將第二大單元中比i大的節(jié)點編號全部減1； （3）修改第二大單元中該節(jié)點所在的單元信息和比該節(jié)點編號要大的那些單元的信息。 重復(fù)以上過程，直至第二大單元中的所有節(jié)點都循環(huán)完畢。,對于三維問題，也可以實現(xiàn)有限元網(wǎng)格的自動分割，只是分塊區(qū)域為六面體或五面體，映射后的區(qū)域為正六面體和直三棱柱。,三 剛度矩陣的帶寬優(yōu)化 為了節(jié)省計算機的存儲空間，需要優(yōu)化總剛矩陣的帶寬。 節(jié)點編號對應(yīng)總剛矩陣中的非零元素的位置。節(jié)點編號變動，總剛矩中非零元素的位置也發(fā)生變動，反之亦然。 帶寬優(yōu)化的原理：通過調(diào)整總剛矩陣中非零元素的位

10、置，對應(yīng)地修改單元信息，從而減少剛陣帶寬。,為此引入鄰接矩陣B，其階數(shù)與剛度矩陣相同，其中的元素非0即1 式中，aij是剛度矩陣中的元素。B也具有帶狀的樣子。 帶寬優(yōu)化有許多實用算法，有些要用到圖論或較深的數(shù)學(xué)知識。下面介紹兩種易于理解的帶寬優(yōu)化方法，它們都是采用變換鄰接矩陣中的行或列的辦法來減少剛陣帶寬的。,1 羅森算法 羅森（Roson）算法的基本思想：找出引起B(yǎng)矩陣最大帶寬的一對頂點，取其中一個與其他頂點互換，看能否減小帶寬。若可以，則進(jìn)行交換。交換后，再繼續(xù)尋找引起當(dāng)時B矩陣最大帶寬的一對頂點，并重復(fù)上述過程。 為此，首先將B中的“頂點”序號紀(jì)錄在一個數(shù)組中，稱為“頂點表”?！绊旤c”是

11、B中每一行離主對角線元素最遠(yuǎn)的那一點。然后在頂點表中： （1）確定B矩陣的帶寬及兩個達(dá)到最大帶寬的第一對頂點。在這對頂點中，若較高編號頂點能與一低編號頂點互換以減小帶寬，則轉(zhuǎn)向（6）；,（2）若較低編號頂點能與一高編號頂點互換以減小帶寬，則轉(zhuǎn)向（6）； （3）若較高編號頂點能與一低編號頂點互換二保持帶寬不變，則轉(zhuǎn)向（6）； （4）若較低編號頂點能與一高編號頂點互換而保持帶寬不變，則轉(zhuǎn)向（6）； （5）如果在（3）、（4）步中已連續(xù)執(zhí)行過規(guī)定的經(jīng)驗次數(shù)，或已交換過的頂點又重新選出來交換，則算法終止；否則，轉(zhuǎn)向（6）； （6）執(zhí)行所指出的頂點交換，即分別交換頂點表的對應(yīng)分量以及B矩陣的行和列，并轉(zhuǎn)

12、向（1）。 可見，整個交換過程都在頂點表中進(jìn)行。,2 阿基茨厄特庫算法 阿基茨(Akyuz)厄特庫(Utku)算法也在1968年提出一種減小平均帶寬的方法，簡稱AU算法。 平均帶寬的定義： AU算法的基本步驟： （1）從B矩陣中取出相鄰的兩行列進(jìn)行交換，并計算平均帶寬。如果滿足下列兩個條件之一，則交換有效：,平均帶寬減少； 平均帶寬保持不變，但有較多元素的行從矩陣中心向外移。 （2）在一個指定的交換循環(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)

13、而平均帶寬不減小，交換運算停止；否則，重復(fù)執(zhí)行(1)、(2)。 一般來講，實用的帶寬優(yōu)化算法并不能使剛度矩陣的帶寬達(dá)到極小化，但并不影響其實用性。 大型的有限元分析軟件包一般都具有帶寬極小化的功能。,四 模型網(wǎng)格的自動生成 目的：檢查所生成的幾何數(shù)據(jù)的正確性。 方法：根據(jù)單元數(shù)組信息（單元和單元節(jié)點關(guān)系）和節(jié)點數(shù)組信息（節(jié)點號和節(jié)點坐標(biāo)關(guān)系），將相鄰節(jié)點連接起來，一個單元一個單元地去形成有限元剖分網(wǎng)格。 注意：所得結(jié)果還要消除隱藏線，否則看不清結(jié)果。這部分內(nèi)容屬于計算機圖形學(xué)，在此不做討論。,8-3 有限元分析的后處理 有限元數(shù)據(jù)后處理的原因： 1）計算結(jié)果是位移和應(yīng)力； 2）數(shù)據(jù)量大，不直觀

14、。 有限元數(shù)據(jù)處理工作包括兩個方面； 1）數(shù)值處理 將有限元結(jié)果轉(zhuǎn)化成工程中常用的形式，或設(shè)計師熟悉的形式； 2 ）圖形處理 用圖形直觀地表示設(shè)計結(jié)果，使設(shè)計結(jié)果一目了然。,一 數(shù)值處理 以應(yīng)力分析為例。有限元位移法計算得到節(jié)點的位移，進(jìn)而計算單元的應(yīng)力： 單元內(nèi)任意一點處的應(yīng)力為 。 上述結(jié)果不論是數(shù)值積分點處的應(yīng)力，節(jié)點處的應(yīng)力，還是單元中心處的應(yīng)力，也不論其精度高低，都不好使用，原因是：我們需要的不是應(yīng)力分量，而是主應(yīng)力或等效應(yīng)力。 說明：單元內(nèi)不同點處的應(yīng)力值的精度是有差別的。一般，數(shù)值積分點處的應(yīng)力值較為準(zhǔn)確。,1 主應(yīng)力計算和等效應(yīng)力計算 工程中往往用一點的主應(yīng)力 來計算或判斷結(jié)構(gòu)

15、的安全程度。對于二維問題，有限元分析可得任意一點的 ，則主應(yīng)力 與它們的關(guān)系為 主應(yīng)力方向,最大剪應(yīng)力和最小剪應(yīng)力為 對于三維應(yīng)力問題，主應(yīng)力可以通過求解下述三次方程得到：,上式的解可表示為 將 按從大到小的順序排列 則 是最大主應(yīng)力， 是最小主應(yīng)力。,從理論上可以證明，三個根都是實根，因此有 否則，一定是有限元計算結(jié)果有誤。 主應(yīng)力的方向可由以下三式中的任意兩式 三式聯(lián)立，得出主應(yīng)力的方向余弦,上式中前兩式與第四式聯(lián)立求解的結(jié)果為 令 ，即可得到 的方向余弦。,另外，要核算所分析的零件強度是否足夠，還要根據(jù)一定的強度理論來計算某點的綜合應(yīng)力。各種強度理論的等效應(yīng)力如下：,2 應(yīng)力修勻 有限元

16、保證位移在單元邊界的連續(xù)性。 但是，應(yīng)力在單元邊界并不連續(xù)。 怎樣處理計算結(jié)果才能使應(yīng)力連續(xù)和直觀呢？ （1）簡單平均或加權(quán)平均 以三角形單元為例。三角形為常應(yīng)力單元，但不同三角形單元的應(yīng)力值不同。 簡單平均： 加權(quán)平均：,（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)達(dá)到最小。 其中,根據(jù) 上式是以修勻后的節(jié)點應(yīng)力向量 為未知量的線性代數(shù)方程組。利用解析法或高斯積分，可以求得 ，進(jìn)而利用單元應(yīng)力插值，可以計算單元內(nèi)任何一點的應(yīng)力值。,單元內(nèi)部

17、分應(yīng)力修勻 只修勻單元中我們感興趣的分量。例如，只修勻分量。這時有,上式為一線性代數(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)為,帶入公式 可以求出未修勻時積分點處的應(yīng)力值 設(shè)修勻后的應(yīng)力在單元內(nèi)按雙線性變化，則修勻應(yīng)力為,令積分點處的應(yīng)力值等于修勻后的應(yīng)力值，則,將上述求得的節(jié)點應(yīng)力，再根據(jù)不同的單元節(jié)點應(yīng)力值予以平均，就可得到最終的節(jié)點應(yīng)力值。 一般情況下，采用局部應(yīng)力磨平處理

18、，可以得到較好的結(jié)果，而計算量是很小的，所以得到廣泛的應(yīng)用。 應(yīng)力全域磨平和子域局部應(yīng)力磨平可參見清華大學(xué)王勖成編著的“有限單元法基本原理和數(shù)值方法”。 二 圖形處理 圖形處理使計算結(jié)果簡潔，直觀，生動。常用以下幾種方法。,1 變形圖 分為兩種： 1）用網(wǎng)格表示內(nèi)部變形（當(dāng)外部變形不大時）； 2 ）用邊界表示外部變形（只關(guān)注外部變形時）。 兩種方法，都把原圖和變形后的圖重疊畫出，以作對比。 繪圖時注意：要放大，否則看不出來。 例如，原節(jié)點坐標(biāo)(x,y)，位移為(u,v)，則繪制變形時該節(jié)點的坐標(biāo)為,2 等值線圖 等值線可以是等溫度線，熱流線，等應(yīng)力線。 方法：追蹤法+插值法。 思路：先找到一個等值點，然后根據(jù)離散節(jié)點坐標(biāo)和節(jié)點的場變量值，求出節(jié)點連線間的等值點位置，并和開始找到的等值點相連。就這樣一直尋找下去，直到邊界。,插值算法：設(shè)兩個節(jié)點i和j的坐標(biāo)為(xi,yi),(xj,yj)，場變量(Ti,Tj)，待求場變量為Ta，則當(dāng) 成立時，在 邊上有一等值點a，其坐標(biāo)為 對于二次插值單元，單元邊界上有三個節(jié)點i,j,k，可采用二次插值或分段線性插值求出等值點。,追蹤算法：從A單元a點到b點，再從b點進(jìn)入B單元，c