第4章 平面問題有限元法-4收斂準(zhǔn)則

1、可以證明，對于一個(gè)給定的位移模式，對于一個(gè)給定的位移模式，其剛度系數(shù)其剛度系數(shù)的數(shù)值比精確值要大。所以，在給定的載荷之下，有限的數(shù)值比精確值要大。所以，在給定的載荷之下，有限元計(jì)算模型的變形將比實(shí)際結(jié)構(gòu)的變形小元計(jì)算模型的變形將比實(shí)際結(jié)構(gòu)的變形小。因而，當(dāng)單元網(wǎng)格分得越來越細(xì)時(shí)，位移的近似解將由下方收斂于位移的近似解將由下方收斂于精確解，即得到精確解，即得到真實(shí)解的下界真實(shí)解的下界。 對于一個(gè)數(shù)值計(jì)算方法，一般總是希望隨著網(wǎng)格的逐步細(xì)分所得到的解答能夠收斂于問題的精確解。根據(jù)前面的分析，我們知道，在有限元分析中，一旦確定了單元的形狀之后，位移模式的選擇將是非常關(guān)鍵的。由于載荷的移置、應(yīng)力矩陣和

2、剛度矩陣的建立等等，都依賴于單元的位移模式，所以，如果所選擇的位移模式與真實(shí)的位移分布有很大的差別，那么就很難獲得良好的數(shù)值解。 位移模式必須能包含單元的常應(yīng)變。位移模式必須能包含單元的常應(yīng)變。每個(gè)單元的應(yīng)變一般都是包含著兩個(gè)部分：一部分是與該單元中各點(diǎn)的坐標(biāo)位置有關(guān)的應(yīng)變（即所謂各點(diǎn)的變應(yīng)變）；另一部分是與位置坐標(biāo)無關(guān)的應(yīng)變（即所謂的常應(yīng)變）。從物理意義上看，為了保證解答的收斂性，要求位移模式必須滿足以下三為了保證解答的收斂性，要求位移模式必須滿足以下三個(gè)條件，即個(gè)條件，即 位移模式必須包含單元的剛體位移。位移模式必須包含單元的剛體位移。也就是說，當(dāng)節(jié)點(diǎn)位移是由某個(gè)剛體位移所引起時(shí)，彈性體內(nèi)

3、將不會產(chǎn)生應(yīng)變。所以，位移模式不但要具有描述單元本身形變的能力，而且還要具有描述由于其它單元形變而通過節(jié)點(diǎn)位移引起單元剛體位移的能力。 例如，三角形三節(jié)點(diǎn)單元位移模式中，常數(shù)項(xiàng)1、4 就是用于提供剛體位移的。 位移模式在單元內(nèi)要連續(xù)、且在相鄰單元之間的位移位移模式在單元內(nèi)要連續(xù)、且在相鄰單元之間的位移必須協(xié)調(diào)。必須協(xié)調(diào)。當(dāng)選擇多項(xiàng)式來構(gòu)成位移模式時(shí)，單元內(nèi)的連續(xù)當(dāng)選擇多項(xiàng)式來構(gòu)成位移模式時(shí)，單元內(nèi)的連續(xù)性要求總是得到滿足的，單元間的位移協(xié)調(diào)性，就是要求單性要求總是得到滿足的，單元間的位移協(xié)調(diào)性，就是要求單元之間既不會出現(xiàn)開裂也不會出現(xiàn)重疊的現(xiàn)象。元之間既不會出現(xiàn)開裂也不會出現(xiàn)重疊的現(xiàn)象。通常，

4、當(dāng)單通常，當(dāng)單元交界面上的位移取決于該交界面上節(jié)點(diǎn)的位移時(shí)，就可以元交界面上的位移取決于該交界面上節(jié)點(diǎn)的位移時(shí)，就可以保證位移的協(xié)調(diào)性。保證位移的協(xié)調(diào)性。 當(dāng)單元尺寸無限縮小時(shí)，每個(gè)單元中的應(yīng)變應(yīng)該趨于常量。當(dāng)單元尺寸無限縮小時(shí)，每個(gè)單元中的應(yīng)變應(yīng)該趨于常量。因此，在位移模式中必須包含有這些常應(yīng)變，否則就不可因此，在位移模式中必須包含有這些常應(yīng)變，否則就不可能使數(shù)值解收斂于正確解。能使數(shù)值解收斂于正確解。 很顯然，三角形三節(jié)點(diǎn)單元位移模式中 ，與2、3、5、6 有關(guān)的線性項(xiàng)就是提供單元中的常應(yīng)變的。 在有限單元法中，把能夠滿足條件在有限單元法中，把能夠滿足條件1 1和和2 2的單元，稱為的單元

5、，稱為完完備單元備單元；滿足條件；滿足條件3 3的單元，叫做的單元，叫做協(xié)調(diào)單元協(xié)調(diào)單元或或保續(xù)單元保續(xù)單元。前。前面討論過的三角形單元和矩形單元，均能同時(shí)滿足上述三個(gè)面討論過的三角形單元和矩形單元，均能同時(shí)滿足上述三個(gè)條件，因此都屬于條件，因此都屬于完備的協(xié)調(diào)單元完備的協(xié)調(diào)單元。在某些梁、板及殼體分。在某些梁、板及殼體分析中，要使單元滿足條件析中，要使單元滿足條件3 3比較困難，所以實(shí)踐中有時(shí)也出比較困難，所以實(shí)踐中有時(shí)也出現(xiàn)一些只滿足條件現(xiàn)一些只滿足條件1 1和和2 2的單元，其收斂性往往也能夠令人滿的單元，其收斂性往往也能夠令人滿意特別是放松條件意特別是放松條件3 3的單元，即的單元，即

6、完備而不協(xié)調(diào)的單元完備而不協(xié)調(diào)的單元，已獲，已獲得了很多成功的應(yīng)用。得了很多成功的應(yīng)用。對于不協(xié)調(diào)單元，其主要的缺點(diǎn)是對于不協(xié)調(diào)單元，其主要的缺點(diǎn)是不不能事先確定其剛度與真實(shí)剛度之間的大小關(guān)系。但值得指出能事先確定其剛度與真實(shí)剛度之間的大小關(guān)系。但值得指出的是，不協(xié)調(diào)單元一般不象協(xié)調(diào)單元那樣剛硬（即比較柔軟的是，不協(xié)調(diào)單元一般不象協(xié)調(diào)單元那樣剛硬（即比較柔軟），因此有可能會比協(xié)調(diào)單元收斂得快。），因此有可能會比協(xié)調(diào)單元收斂得快。在選擇多項(xiàng)式作為單元的位移模式時(shí)，其階次的確定，在選擇多項(xiàng)式作為單元的位移模式時(shí)，其階次的確定，要考慮解答的收斂性要考慮解答的收斂性，即單元的完備性和協(xié)調(diào)性要求。實(shí)踐，

7、即單元的完備性和協(xié)調(diào)性要求。實(shí)踐證明，這兩項(xiàng)確實(shí)是所要考慮的重要因素，但并不是唯一的證明，這兩項(xiàng)確實(shí)是所要考慮的重要因素，但并不是唯一的因素。選擇多項(xiàng)式位移模式階次時(shí)，需要考慮的另一個(gè)因素因素。選擇多項(xiàng)式位移模式階次時(shí)，需要考慮的另一個(gè)因素是，是，所選的模式應(yīng)該與局部坐標(biāo)系的方位無關(guān)所選的模式應(yīng)該與局部坐標(biāo)系的方位無關(guān)，這一性質(zhì)稱，這一性質(zhì)稱為為幾何各向同性幾何各向同性。對于線性多項(xiàng)式，各向同性的要求通常就。對于線性多項(xiàng)式，各向同性的要求通常就等價(jià)于位移模式必須包含常應(yīng)變狀態(tài)。對于高次位移模式，等價(jià)于位移模式必須包含常應(yīng)變狀態(tài)。對于高次位移模式，就是不應(yīng)該有一個(gè)偏惠的坐標(biāo)方向，也就是位移形式不

8、應(yīng)該就是不應(yīng)該有一個(gè)偏惠的坐標(biāo)方向，也就是位移形式不應(yīng)該隨局部坐標(biāo)的更換而改變。經(jīng)驗(yàn)證明，實(shí)現(xiàn)幾何各向同性的隨局部坐標(biāo)的更換而改變。經(jīng)驗(yàn)證明，實(shí)現(xiàn)幾何各向同性的一種有效方法是，根據(jù)如圖所示的一種有效方法是，根據(jù)如圖所示的巴斯卡三角形巴斯卡三角形來選擇二維來選擇二維多項(xiàng)式的各項(xiàng)。在二維多項(xiàng)式中，如果包含有對稱軸一邊的多項(xiàng)式的各項(xiàng)。在二維多項(xiàng)式中，如果包含有對稱軸一邊的某一項(xiàng)，那么就必須同時(shí)包含有另一邊的對稱項(xiàng)。某一項(xiàng)，那么就必須同時(shí)包含有另一邊的對稱項(xiàng)。選擇多項(xiàng)式位移模式時(shí)，還應(yīng)該要考慮的一個(gè)因素：選擇多項(xiàng)式位移模式時(shí)，還應(yīng)該要考慮的一個(gè)因素：多多項(xiàng)式中的項(xiàng)數(shù)必須等于或稍大于單元邊界上的外節(jié)點(diǎn)的

9、自由項(xiàng)式中的項(xiàng)數(shù)必須等于或稍大于單元邊界上的外節(jié)點(diǎn)的自由度數(shù)度數(shù)。通常是取項(xiàng)數(shù)與單元的外節(jié)點(diǎn)的自由度數(shù)相等，取過。通常是取項(xiàng)數(shù)與單元的外節(jié)點(diǎn)的自由度數(shù)相等，取過多的項(xiàng)數(shù)是不恰當(dāng)?shù)?。多的?xiàng)數(shù)是不恰當(dāng)?shù)摹?22322343223454322345xyxxyyxx yxyyxx yx yxyyxx yx yx yxyy常數(shù)項(xiàng)線性項(xiàng)二次項(xiàng)三次項(xiàng)四次項(xiàng)五次項(xiàng)對稱軸巴斯卡三角形根據(jù)前面的討論，現(xiàn)以三角形常應(yīng)變單元為例來說明應(yīng)用有限元法求解彈性力學(xué)平面問題的具體步驟：力學(xué)模型的確定根據(jù)工程實(shí)際情況確定問題的力學(xué)模型，并按一定比例繪制結(jié)構(gòu)圖、注明尺寸、載荷和約束情況等。將計(jì)算對象進(jìn)行離散化，即彈性體劃分為許多

10、三角形單元，并對節(jié)點(diǎn)進(jìn)行編號。確定全部節(jié)點(diǎn)的坐標(biāo)值，對單元進(jìn)行編號，并列出各單元三個(gè)節(jié)點(diǎn)的節(jié)點(diǎn)號。 計(jì)算載荷的等效節(jié)點(diǎn)力（要求的輸入信息）。 由各單元的常數(shù)bi 、ci 、bj 、cj 、bm 、cm 及行列式2，計(jì)算單元剛度矩陣。 組集整體剛度矩陣，即形成總剛的非零子矩陣。 處理約束，消除剛體位移。 求解線性方程組，得到節(jié)點(diǎn)位移。 計(jì)算應(yīng)力矩陣，求得單元應(yīng)力，并根據(jù)需要計(jì)算主應(yīng)力和主方向。 整理計(jì)算結(jié)果（后處理部分）。 為了提高有限元分析計(jì)算的效率、達(dá)到一定的精度，應(yīng)該注意以下幾個(gè)方面。 一一. . 對稱性的利用對稱性的利用 在劃分單元之前，有必要先研究一下計(jì)算對象的對稱或反對稱的情況，以便

11、確定是取整個(gè)物體，還是部分物體作為計(jì)算模型。 例如，圖4-11(a)所示受純彎曲的梁，其結(jié)構(gòu)對于x、y軸都是幾何對稱的，而所受的載荷則是對于x軸對稱，對于x軸反對稱。可知，梁的應(yīng)力和變形也將具有同樣的對稱特性，所以只需取1/4梁進(jìn)行計(jì)算即可。取分離體如圖4-11(b)所示，對于其它部分結(jié)構(gòu)對此分離體的影響，可以作相應(yīng)的處理，即對處于y軸對稱面內(nèi)各節(jié)點(diǎn)的x方向位移都設(shè)置為零，而對于在x軸反對稱面上的各節(jié)點(diǎn)的x方向位移也都設(shè)置為零。這些條件就等價(jià)于在圖4-11(b)中相應(yīng)節(jié)點(diǎn)位置處施加約束，圖中o點(diǎn)y方向施加的約束是為了消除剛體位移。RRRRyxooyxR(a)(b)節(jié)點(diǎn)的布置是與單元的劃分互相聯(lián)

12、系的。通常，集中載荷的作用點(diǎn)、分布載荷強(qiáng)度的突變點(diǎn)，分布載荷與自由邊界的分界點(diǎn)、支承點(diǎn)等都應(yīng)該取為節(jié)點(diǎn)。并且，當(dāng)物體是由不二二. . 節(jié)點(diǎn)的選擇及單元的劃分節(jié)點(diǎn)的選擇及單元的劃分圖4-11節(jié)點(diǎn)的多少及其分布的疏密程度（即單元的大?。?，一般要根據(jù)所要求的計(jì)算精度等方面來綜合考慮。從計(jì)算結(jié)果的精度上講，當(dāng)然是單元越小越好，但計(jì)算所需要的時(shí)間也要大大增加。另外，在微機(jī)上進(jìn)行有限元分析時(shí)，還要考慮計(jì)算機(jī)的容量。因此，在保證計(jì)算精度的前提下，應(yīng)力求采用較少的單元。為了減少單元(a) (b)圖4-12同的材料組成時(shí)，厚度不同或材料不同的部分，也應(yīng)該劃分為不同的單元。在進(jìn)行節(jié)點(diǎn)編號時(shí)，應(yīng)該注意要盡量使同一單

13、元的相鄰節(jié)點(diǎn)的號碼差盡可能地小，以便最大限度地縮小剛度矩陣的帶寬，節(jié)省存儲、提高計(jì)算效率。如前所述，平面問題的半帶寬為B =2(d+1)，在劃分單元時(shí)，對于應(yīng)力變化梯度較大的部位單元可小一些，而在應(yīng)力變化比較平緩的區(qū)域可以劃分得粗一些。還有一點(diǎn)值得注意的是，單元各邊的長度不要相差太大，以免出現(xiàn)過大的計(jì)算誤差或出現(xiàn)病態(tài)矩陣。例如，圖4-12所示的(a)、(b)兩種單元劃分，雖然都是同樣的四個(gè)節(jié)點(diǎn)，但(a)的劃分方式顯然要比(b)的方式好。三三. . 節(jié)點(diǎn)的編號節(jié)點(diǎn)的編號1 2 3 4 5 6 71 3 5 7 9 11 138 9 10 11 12 13 142 4 6 8 10 12 14若采

14、取帶寬壓縮存儲，則整體剛度矩陣的存儲量N 最多為N =2nB = 4n(d+1)其中：d為相鄰節(jié)點(diǎn)的最大差值，n為節(jié)點(diǎn)總數(shù)。例如在圖4-13中，(a)與(b)的單元劃分相同，且節(jié)點(diǎn)總數(shù)都等于14，但兩者的節(jié)點(diǎn)編號方式卻完全不同。(a)是按長邊進(jìn)行編號，d =7，N =488；而(b)是按短邊進(jìn)行編號，d =2，N =168。顯然(b)的編號方式可比(a)的編號方式節(jié)省280個(gè)存儲單元。(a) (b)圖4-13四四. . 單元節(jié)點(diǎn)單元節(jié)點(diǎn)i i、j j、m m的次序的次序在前面章節(jié)中，我們曾指出，為了在計(jì)算中保證單元的面積 不會出現(xiàn)負(fù)值，節(jié)點(diǎn)i、j、m的編號次序必須是逆時(shí)針方向。事實(shí)上，節(jié)點(diǎn)i、

15、j、m的編號次序是可以任意安排的，只要在計(jì)算剛度矩陣的各元素時(shí)，對取絕對值，即可得到正確的計(jì)算結(jié)果。在實(shí)際計(jì)算時(shí)，應(yīng)該注意所選有限元分析軟件的使用要求。五五. . 邊界條件的處理及整體剛度矩陣的修正邊界條件的處理及整體剛度矩陣的修正 在前面討論整體剛度矩陣時(shí)，已經(jīng)提到，整體剛度矩陣的奇異性可以提高考慮邊界約束條件來排除彈性體的剛體位移，以達(dá)到求解的目的。一般情況下，求解的問題，其邊界往往已有一點(diǎn)的位移約束條件，本身已排除了剛體運(yùn)動的可能性。否則的話，就必須適當(dāng)指定某些節(jié)點(diǎn)的位移值，以避免出現(xiàn)剛體位移。這里介紹兩種比較簡單的引入已知節(jié)點(diǎn)位移的方法，這兩種方法都可保持原K矩陣的稀疏、帶狀和對稱等特

16、性。 下面我們來實(shí)際考察一個(gè)只有四個(gè)方程的簡單例子。 保持方程組為2n2n系統(tǒng)，僅對K和R進(jìn)行修正。例如，若指定節(jié)點(diǎn)i在方向y的位移為vi ，則令K中的元素k2i, 2i 為1，而第2i行和第2i列的其余元素都為零。R中的第2i個(gè)元素則用位移vi 的已知值代入，R中的其它各行元素均減去已知節(jié)點(diǎn)位移的指定值和原來K中該行的相應(yīng)列元素的乘積。KKKKKKKKKKKKKKKKuvuvRRRR11121314212223243132333441424344112212341000000010002224424411221221123334411433KKKKuvuvRKKRKK假定該系統(tǒng)中節(jié)點(diǎn)位移u1

17、 和u2分別被指定為當(dāng)引入這些節(jié)點(diǎn)的已知位移之后，方程(a)就變成 然后，就用這組維數(shù)不變的方程來求解所有的節(jié)點(diǎn)位移。顯然，其解答仍為原方程(a)的解答。 u1 = 1 ， u2 = 2KKKKKKKKKKKKKKKKuvuvKRKR111512131421222324313233153441424344112211115233315410101010KuK vK uK vK11151121132142111151010將K中與指定的節(jié)點(diǎn)位移有關(guān)的主對角元素乘上一個(gè)大數(shù)，如1015，同時(shí)將R中的對應(yīng)元素?fù)Q成指定的節(jié)點(diǎn)位移值與該大數(shù)的乘積。實(shí)際上，這種方法就是使K中相應(yīng)行的修正項(xiàng)遠(yuǎn)大于非修正項(xiàng)。

18、 若把此方法用于上面的例子，則方程(a)就變成 事實(shí)上，該方程組的第一個(gè)方程為 圖4-14所示為一厚度t=1cm的均質(zhì)正方形薄板，上下受均勻拉力q=106N/m，材料彈性模量為E，泊松比 ，不記自重，試用有限元法求其應(yīng)力分量。3/1123421x圖 4-15y2q2q2myxq=106N/m圖 4-14011213132321yybyybyyb解：解：.力學(xué)模型的確定.結(jié)構(gòu)離散由于此結(jié)構(gòu)長、寬遠(yuǎn)大于厚度，而載荷作用于板平面內(nèi)，且沿板厚均勻分布，故可按平面應(yīng)力問題處理，考慮到結(jié)構(gòu)和載荷的對稱性，可取結(jié)構(gòu)的1/4來研究。該1/4結(jié)構(gòu)被離散為兩個(gè)三角形單元，節(jié)點(diǎn)編號，單元劃分及取坐標(biāo)如圖4-15所示

19、，其各節(jié)點(diǎn)的坐標(biāo)值見表4-1。.求單元的剛度矩陣計(jì)算單元的節(jié)點(diǎn)坐標(biāo)差及單元面積單元（i、j、m1，2，3）節(jié)點(diǎn)坐標(biāo)xy表4-1 計(jì)算各單元的剛度矩陣 先計(jì)算用到的常數(shù)211011212111023321213132321cbcbxxcxxcxxc89)1 (2169)1 (4312122EEEEt 310011691131000131103110310131003111169111EEK代入可得:34323234169;031310169;3131311169122113112EKEKEK10031169;1313131169133123EKEK所以單元1的剛度矩陣為: 10311313431

20、31323403131313131031101169133132131123122121113112111166稱對EKKKKKKKKKK123123由于單元2若按3、4、1對應(yīng)單元1的1、2、3排碼時(shí),則這兩個(gè)單元剛度矩陣內(nèi)容完全一樣，故有: 1031131343131323403131313131031101169266稱對EK341341組集整體剛度矩陣 由于Krs=KsrT，又單元1和單元2的節(jié)點(diǎn)號按1、2、3對應(yīng)3、4、1，則可得： 24424324121331322131122121211188KKKKKKKKKK0110163111316310031631132131311122

21、43121233111EKKKEKKKEKKTT按剛度集成法可得整體剛度矩陣為：31111630110163300116331111634224163214123241213113231211133123214132244122EKKKEKKKEKKEKKKEKKTTT所以組集的整體剛度矩陣為： 42110031413001140310240120421141340416388稱對EK先求出各單元的應(yīng)力矩陣S1、S2，然后再求得各單元的應(yīng)力分量：6. 計(jì)算各單元應(yīng)力矩陣，求出各單元應(yīng)力 01003/8083/3/03/0001111030310110130383111qqEqEqEqESxyy

22、x 0100)3/(808300/0/3/01111030310110130383222qEqEEqEqEqESxyyx單元應(yīng)力可看作是單元形心處的應(yīng)力值。7. 引入約束條件，修改剛度方程并求解根據(jù)約束條件：u1 =v1=0；v2=0；u4=0和等效節(jié)點(diǎn)力列陣： ， 并 代 入 剛 度 方 程 ： ，劃去K中與0位移相對應(yīng)的1，2，4，7的行和列，則剛度方程變?yōu)椋?TqqF2/02/00000 FK2/2/0041104014141634332qqvvuuETTEqEqEqEqvvuu/3/3/4332 TEq1013/103/100/求解上面方程組可得出節(jié)點(diǎn)位移為：所以 圖4-16所示為一平

23、面應(yīng)力問題離散化以后的結(jié)構(gòu)圖，其中圖（a）為離散化后的總體結(jié)構(gòu)，圖（b）為單元1，2，3，4的結(jié)構(gòu)，圖（c）為單元3的結(jié)構(gòu)。用有限元法計(jì)算節(jié)點(diǎn)位移、單元應(yīng)變及單元應(yīng)力（為簡便起見，取泊松比 ，單元厚度t=1）。0 xyp1234651234a3ijmivaiXiujXjuiYjvjYmvmYmXmuaa1,2,4ijmiXiviujXjuiYjvjYmvmYmXmu 圖 4-16 計(jì)算實(shí)例2的結(jié)構(gòu)圖首先求確定各單元剛度所需的系數(shù) 及面積A，對于單元1，2，4有：mjimjicccbbb,2/0, 02aAcacacababbmjimji解：解：對于單元3有:2/, 0, 0,2aAacacca

24、bbabmjimji 其次，求出各單元的單元剛度矩陣。對于1，2，4單元，其單元剛度矩陣為：10110102020010312112130100202010110144, 2, 1Ekijmijm 各單元的節(jié)點(diǎn)編號與總體結(jié)構(gòu)的總編號之間的對應(yīng)關(guān)系見表4-2。 31211013011220200011011011011002000243Ek 對于單元3，其單元剛度矩陣為：ijmijm 各單元節(jié)點(diǎn)號與總體節(jié)點(diǎn)號對應(yīng)表各單元節(jié)點(diǎn)號與總體節(jié)點(diǎn)號對應(yīng)表 單元號 1 2 3 4 節(jié)點(diǎn)號 節(jié) 點(diǎn) 總 編 號 I 1 2 2 3 j 2 4 5 5 m 3 5 3 6表表4-2 將各單元剛度矩陣按節(jié)點(diǎn)總數(shù)及相

25、應(yīng)的節(jié)點(diǎn)號關(guān)系擴(kuò)充成12*12矩陣,分別如下: 654321000000000000000000000000000000000000000000000000000000000000000000000000000000101101000000020200000000103121000000121301000000002020000000101101465432111212Ek6543210000000000000000000000000010110001000002020000000010310021000012130001000000000000000000000000000000200020

26、00001011000100000000000000000000000000465432121212Ek654321000000000000000000000000002000200100000100111000000000000000000000000000002100311000000100131200000100111000000000020200000000000000000000000000465432131212Ek654321101100010000020200000000103100210000121300010000000000000000000000000000002000200000101100010000000000000000000000000000000000000000000000000000465432141212Ek將擴(kuò)充后的各單元剛度矩陣相加,得總體剛度矩陣K，即： 65432110110001000002020000000010611141010012160212100000103100210000121300010000410061210110120016140000012021612100101