有限元分析基礎(chǔ)平面問題的有限元

第三章平面問題的有限元§3.1三角形單元§3.2矩形單元§3.3平面等參元§3.4平面單元應(yīng)用比較分析第三章平面問題的有限元3.1三角形單元§3.1.1位移函數(shù)§3.1.2應(yīng)變與應(yīng)力矩陣§3.1.3形函數(shù)的性質(zhì)§3.1.4單元剛度矩陣及其特點§3.1.5等效節(jié)點載荷§3.1.6整體剛度矩陣及其特點§3.1.7邊界條件§3.1.8算例§3.1.9收斂準(zhǔn)則§3.1.10從能量原理推導(dǎo)剛度矩陣本節(jié)主要研究三節(jié)點三角形單元，首先建立以單元節(jié)點位移表示單元內(nèi)各點位移的關(guān)系式。設(shè)單元e的節(jié)點編號為i、j、m，如圖所示。由彈性力學(xué)平面問題可知，每個節(jié)點在其單元平面內(nèi)的位移有兩個分量，所以整個三角形單元將有六個節(jié)點位移分量，即六個自由度。用列陣可表示為式中的子向量為

（i、j、m輪換）式中：——節(jié)點i在x軸方向的位移分量；

——節(jié)點i在y軸方向的位移分量。3.1.1位移函數(shù)在有限單元法中，雖然是用離散化模型來代替原來的連續(xù)體，但每一個單元體仍是一個彈性體，所以在其內(nèi)部依然是符合彈性力學(xué)基本假設(shè)的，彈性力學(xué)的基本方程在每個單元內(nèi)部同樣適用。從彈性力學(xué)平面問題的解析解法中可知，如果彈性體內(nèi)的位移分量函數(shù)已知，則應(yīng)變分量和應(yīng)力分量也就確定了。但是，如果只知道彈性體中某幾個點的位移分量的值，那么就不能直接求得應(yīng)變分量和應(yīng)力分量。因此，在進行有限元分析時，必須先假定一個位移模式。由于在彈性體內(nèi)，各點的位移變化情況非常復(fù)雜，很難在整個彈性體內(nèi)選取一個恰當(dāng)?shù)奈灰坪瘮?shù)來表示位移的復(fù)雜變化，但是如果將整個區(qū)域分割成許多小單元，那么在每個單元的局部范圍內(nèi)就可以采用比較簡單的函數(shù)來近似地表示單元的真實位移，將各單元的位移模式連接起來，便可近似地表示整個區(qū)域的真實位移函數(shù)。這種化繁為簡、聯(lián)合局部逼近整體的思想，正是有限單元法的絕妙之處。3.1.1位移函數(shù)基于上述思想，可以選擇一個單元位移模式，單元內(nèi)各點的位移可按此位移模式由單元節(jié)點位移通過插值而獲得。假設(shè)三節(jié)點三角形單元內(nèi)任意點的位移是x、y的線性函數(shù)，如下是待定常數(shù)，因三角形單元共有六個自由度，且位移函數(shù)在三個節(jié)點處的數(shù)值應(yīng)該等于這些點處的位移分量的數(shù)值。假設(shè)節(jié)點i、j、m的坐標(biāo)分別為（式中～、）、（）、（），代入位移函數(shù)，得3.1.1位移函數(shù)從上式中可以解得式中：——三角形i、j、m的面積。為保證求得的面積為正值，節(jié)點i、j、m的編排次序必須是逆時針方向，如圖所示。的計算公式，如下3.1.1位移函數(shù)將式所求得得系數(shù)代入位移函數(shù)，得式中：（i、j、m輪換）3.1.1位移函數(shù)同理可以求出單元內(nèi)任意點的y方向位移，如下令

（i、j、m輪換）則單元內(nèi)任意點的位移可以表示成如下形式式中：——單元內(nèi)位移矢量3.1.1位移函數(shù)位移函數(shù)可以進一步寫成如下形式式中：

（i、j、m輪換）是坐標(biāo)的函數(shù)，它們反映了單元的位移狀態(tài)，所以一般稱之為形狀函數(shù)，簡稱形函數(shù)。顯然，形函數(shù)決定了單元內(nèi)的“位移模式”，反映了i節(jié)點位移對單元內(nèi)任意點位移的貢獻。3.1.1位移函數(shù)有了單元的位移模式，就可以利用平面問題的幾何方程求出單元內(nèi)任意點的應(yīng)變，式中：——單元應(yīng)變列陣；

——單元應(yīng)變轉(zhuǎn)換矩陣，或稱單元幾何矩陣。3.1.2應(yīng)變與應(yīng)力矩陣單元應(yīng)變轉(zhuǎn)換矩陣，為3×6矩陣，其子矩陣為（i、j、m輪換）可以看出，表示i節(jié)點位移對單元應(yīng)變的貢獻率。當(dāng)單元確定后，也就確定了，此時單元內(nèi)的應(yīng)變僅依賴于節(jié)點的位移。對于三節(jié)點三角形單元，由于面積和

等都是常量，所以單元應(yīng)變轉(zhuǎn)換矩陣中的諸元素都是常量，因而單元中各點的應(yīng)變分量也都是常量，即單元內(nèi)各點的應(yīng)變相同，通常稱這種單元為常應(yīng)變單元。、3.1.2應(yīng)變與應(yīng)力矩陣求得應(yīng)變之后，再將應(yīng)變代入物理方程，即公式（2.19），便可推導(dǎo)出以節(jié)點位移表示的應(yīng)力，如下式中：——單元應(yīng)力矩陣。令則式中：——單元應(yīng)力轉(zhuǎn)換矩陣。3.1.2應(yīng)變與應(yīng)力矩陣的子矩陣為單元應(yīng)力轉(zhuǎn)換矩陣可以看出，表示i節(jié)點位移對單元應(yīng)力的貢獻率。當(dāng)單元確定后，也就確定了，此時單元內(nèi)的應(yīng)力僅依賴于節(jié)點的位移。對于三節(jié)點三角形單元，由于面積和

等都是常量，所以單元應(yīng)力轉(zhuǎn)換矩陣中的諸元素都是常量，因而單元中各點的應(yīng)力分量也都是常量，即單元內(nèi)各點的應(yīng)力相同，通常稱這種單元為常應(yīng)力單元。、3.1.2應(yīng)變與應(yīng)力矩陣可見，對于三節(jié)點三角形單元，由于所選取的位移模式是線性的，因而導(dǎo)致單元內(nèi)各點的應(yīng)變、應(yīng)力相同，但相鄰單元將具有不同的應(yīng)力和應(yīng)變，即在單元的公共邊界上應(yīng)力和應(yīng)變的值將會有突變，但位移卻是連續(xù)的。3.1.2應(yīng)變與應(yīng)力矩陣由前文可知，形函數(shù)決定了單元內(nèi)的“位移模式”，反映了i節(jié)點位移對單元內(nèi)任意點位移的貢獻，形函數(shù)具有如下性質(zhì)：1.形函數(shù)在各單元節(jié)點上的值，具有“本點是1、它點為零”的性質(zhì)即

（i、j、m輪換）證明如下：3.1.3形函數(shù)的性質(zhì)在節(jié)點j、m上類似有在節(jié)點i上3.1.3形函數(shù)的性質(zhì)2.在單元內(nèi)任意點上，三個形函數(shù)之和等于1即證明：將單元內(nèi)任意點的坐標(biāo)代入上式，得，因此容易得到此結(jié)論。由此可見，三個形函數(shù)中只有二個是獨立的。可以證明，，3.1.3形函數(shù)的性質(zhì)例如：假設(shè)單元各節(jié)點的位移均相同并等于，則根據(jù)此性質(zhì)可得單元內(nèi)任意點的位移為，顯然，該性質(zhì)反映了單元剛體位移，即當(dāng)單元做剛體運動時，單元內(nèi)任意點的位移均等于剛體位移。3.1.3形函數(shù)的性質(zhì)3.三角形單元任意一條邊上的形函數(shù)，僅與該邊的兩端節(jié)點坐標(biāo)有關(guān)。例如在ij邊上，有事實上，因ij

邊的直線方程方程為3.1.3形函數(shù)的性質(zhì)代入及，得則3.1.3形函數(shù)的性質(zhì)另外，由可以求得利用形函數(shù)的這一性質(zhì)可以證明，相鄰單元的位移分別進行線性插值之后，在其公共邊上將是連續(xù)的。例如，對圖3.2所示的單元具有ij公共邊，可知，在ij邊上有這樣，不論按哪個單元來計算，公共邊ij上的位移均由下式表示3.1.3形函數(shù)的性質(zhì)由此可見，在公共邊上的位移u、v將完全由公共邊上的兩個節(jié)點i、j的位移所確定，因而相鄰單元的位移是連續(xù)的。例3.1求如圖3.3所示等邊三角形的形函數(shù)。三角形的面積將以上數(shù)據(jù)代入式（3.10），得3.1.3形函數(shù)的性質(zhì)①將各節(jié)點坐標(biāo)代入以上三式，有即滿足“本點是1、它點為零”的性質(zhì)。3.1.3形函數(shù)的性質(zhì)即滿足“在單元內(nèi)任意點上，三個形函數(shù)之和等于1”的性質(zhì)。②將三個形函數(shù)相加，得，即在ij邊上任意位置的位移，將完全由節(jié)點i、j的位移所確定。取其它邊同樣可以得出其結(jié)論，即滿足“三角形單元任意一條邊上的形函數(shù)，僅與該邊的兩端節(jié)點坐標(biāo)有關(guān)”的性質(zhì)。③以ij邊為例，在ij邊上，則3.1.3形函數(shù)的性質(zhì)單元在節(jié)點處受力，單元會發(fā)生變形，也就是說單元在節(jié)點處受到的力與單元節(jié)點位移之間有必然的聯(lián)系。單元間正是通過與節(jié)點間的相互作用力連接起來成為整體，而一個單元仍是一個彈性體，如果將每一個單元的受力－位移關(guān)系找到，則整體的受力－位移關(guān)系也容易清楚。為了推導(dǎo)單元的節(jié)點力和節(jié)點位移之間的關(guān)系，可應(yīng)用虛位移原理對單元進行分析。單元在節(jié)點處受到的力稱為單元節(jié)點力，單元在節(jié)點力的作用下處于平衡，節(jié)點力可以表示為

式中：——節(jié)點i沿x方向的節(jié)點力分量；

——節(jié)點i沿y方向的節(jié)點力分量。3.1.4單元剛度矩陣及其特點單元節(jié)點的虛位移為相應(yīng)的單元內(nèi)的虛位移場為

單元內(nèi)的虛應(yīng)變?yōu)?/p>

設(shè)單元僅在節(jié)點上受力，則單元節(jié)點力在節(jié)點虛位移上的虛功為

3.1.4單元剛度矩陣及其特點單元內(nèi)的應(yīng)力在虛應(yīng)變上所做的功為式中：——在體積上積分。令單元剛度矩陣為由虛功原理可知，式（3.27）和（3.28）相等，再由虛位移的任意性，可得3.1.4單元剛度矩陣及其特點上式為單元剛度矩陣的定義式，其推導(dǎo)過程與形函數(shù)的具體表達形式、節(jié)點個數(shù)均無關(guān)，這說明該表達式具有普遍意義，則上式是單元節(jié)點力與單元節(jié)點位移之間的關(guān)系式，即單元節(jié)點平衡方程，或稱單元剛度方程，具有普遍意義。對于三節(jié)點三角形單元，由于、中的元素都是常量，這里假設(shè)單元厚度均勻，則式中：——單元面積；

——單元厚度。3.1.4單元剛度矩陣及其特點上式進一步可以表示為將幾何矩陣式（3.16）和彈性矩陣式（2.19）代入，可得每個子塊的矩陣表達式，如下3.1.4單元剛度矩陣及其特點單元剛度矩陣具有如下特點：1.中的每一個元素都是一個剛度系數(shù)假設(shè)單元的節(jié)點位移為，由，得到節(jié)點力為

3.1.4單元剛度矩陣及其特點由此可得：表示i節(jié)點在水平方向產(chǎn)生單位位移時，在節(jié)點i的水平方向上需要施加的節(jié)點力；表示i節(jié)點在水平方向產(chǎn)生單位位移時，在節(jié)點i的垂直方向上需要施加的節(jié)點力。選擇不同的單元節(jié)點位移，可以得到單元剛度矩陣中每個元素的物理含義，如下：表示s節(jié)點在水平方向產(chǎn)生單位位移時，在節(jié)點r的水平方向上需要施加的節(jié)點。表示s節(jié)點在水平方向產(chǎn)生單位位移時，在節(jié)點r的垂直方向上需要施加的節(jié)點力；表示s節(jié)點在垂直方向產(chǎn)生單位位移時，在節(jié)點r的水平方向上需要施加的節(jié)點力；表示s節(jié)點在水平方向產(chǎn)生單位位移時，在節(jié)點r的垂直方向上需要施加的節(jié)點力；因此，中的每一個元素都是一個剛度系數(shù)，表示單位節(jié)點位移分量所引起的節(jié)點力分量，中的每一個子塊表示單元某個節(jié)點位移矢量對單元某個節(jié)點力矢量的貢獻率。3.1.4單元剛度矩陣及其特點2.是對稱矩陣由式（3.32）可知，是對稱矩陣。根據(jù)對稱性，可以減少計算量和存儲量。

3.是奇異矩陣，即根據(jù)這一性質(zhì)，當(dāng)已知單元節(jié)點位移時，可以從式（3.31）中求出單元節(jié)點力。反之，由于單元剛陣奇異，不存在逆矩陣，因此，當(dāng)已知單元節(jié)點力時不能求出單元上的節(jié)點位移。從物體變形的實際情況來說，單元剛陣奇異是必須的，因為單元除了產(chǎn)生變形外，還會產(chǎn)生剛體位移，即僅僅依靠節(jié)點力是無法唯一地確定剛體位移的。事實上，當(dāng)單元的節(jié)點力為零時，單元仍可做剛體運動。3.1.4單元剛度矩陣及其特點例如，假定單元產(chǎn)生了x方向的單位剛體位即并假設(shè)此時對應(yīng)的單元節(jié)點力為零，則由

得可以得到，在單元剛度矩陣中1，3，5列中對應(yīng)行的系數(shù)相加為零，由行列式的性質(zhì)可知，。3.1.4單元剛度矩陣及其特點

4.單元的剛度不隨單元或坐標(biāo)軸的平行移動而改變單元的剛度取決于單元的形狀、大小、方向和彈性系數(shù)，而與單元的位置無關(guān)，即不隨單元或坐標(biāo)軸的平行移動而改變。因此，只要單元的形狀、大小、方向和彈性系數(shù)相同，無論單元出現(xiàn)在任何位置均有相同的單元剛度矩陣，根據(jù)對稱性，可以減少計算量。，厚度為t。例3.2求圖3.4所示等腰直角三角形單元的剛度矩陣，設(shè)解根據(jù)1（0，a）、2（0，0）、3（a，0）的坐標(biāo)和式（3.8）可得3.1.4單元剛度矩陣及其特點單元面積，將以上各數(shù)據(jù)代入幾何矩陣式（3.16）為由于，則平面應(yīng)力與平面應(yīng)變問題的彈性矩陣相同，為3.1.4單元剛度矩陣及其特點由式（3.18）可得應(yīng)力轉(zhuǎn)換矩陣，如下則單元剛度矩陣為3.1.4單元剛度矩陣及其特點從以上求得的結(jié)果可以看出，單元剛度矩陣是對稱的。1、3、5列相加等于零，2、4、6列相加也等于零，單元剛度矩陣是奇異的。另外在幾何矩陣中所涉及到的元素均為“坐標(biāo)差量”，而與坐標(biāo)的具體值無關(guān)，從而也說明了單元剛度矩陣不隨單元或坐標(biāo)軸的平行移動而改變。1）獲取單元節(jié)點信息；2）計算單元面積及參數(shù)、等；3）計算幾何矩陣中各元素；4）根據(jù)問題類型（平面應(yīng)力、平面應(yīng)變）計算彈性矩陣；5）計算應(yīng)力轉(zhuǎn)換矩陣；6）計算單元剛度矩陣。從單元剛度矩陣的計算公式可以看出，在求解單元剛度矩陣時，可以不必求得應(yīng)力轉(zhuǎn)換矩陣，但在后期求解單元應(yīng)力時要用到該矩陣，所以在此處一并求出，也可以根據(jù)實際情況具體考慮。3.1.4單元剛度矩陣及其特點實際上作用在物體上的外力可以直接作用在節(jié)點上，也可以不作用在節(jié)點上，然而在有限元法中要求作用在物體上的各種外力必須用作用在節(jié)點上的力表示。這一點體現(xiàn)了有限法中“離散化”這一概念，即將連續(xù)體離散成只有在節(jié)點處相連的單元體，單元與單元之間的聯(lián)系只能通過節(jié)點。如前文所述，單元內(nèi)任意點的位移、應(yīng)力、應(yīng)變等變量最終都用單元節(jié)點位移來表示。同樣，作用在物體上的各種外力也必須用作用在節(jié)點上的力表示，這一過程稱為外力的靜力等效移置，所得到的節(jié)點力稱為等效節(jié)點力。

3.1.5等效節(jié)點載荷靜力等效原則：

對于剛體來說，所謂靜力等效原則就是單元上原有的外力系和將外力系向各節(jié)點移置所得的等效節(jié)點力，二者向同一點簡化應(yīng)具有相同的主失和主矩；對于彈性體來說，所謂靜力等效原則就是指單元上的外力系和將該力系向各節(jié)點移置后的等效節(jié)點力，二者在虛位移上的虛功相等，也即外力作用在單元上所引起的變形能和移置后的等效節(jié)點力在單元上引起的變形能相等，在一定的位移模式下這種移置是唯一的。3.1.5等效節(jié)點載荷1.彈性體靜力等效原則——虛功原理虛功等效的思路是：就一個單元來說，把作用在單元上的外力移置到節(jié)點上后，應(yīng)當(dāng)與原來的實際外力所作虛功等效。其計算方法是：對任意允許的微小虛位移，原外力所作虛功等于移置后的等效節(jié)點力所作虛功。即式中：——節(jié)點虛位移；——單元虛位移場；——集中力作用處產(chǎn)生的虛位移。3.1.5等效節(jié)點載荷等式左邊表示等效節(jié)點力在節(jié)點虛位移上所作虛功；右邊第一項表示作用在單元任意點c處的集中力所作虛功、第二項表示表面力所作的虛功、第三項表示體積力所作的虛功。將代入上式，并考慮的任意性，得式中：——由集中力等效的單元等效節(jié)點力；——由表面力等效的單元等效節(jié)點力；——由體積力等效的單元等效節(jié)點力；——形函數(shù)在集中力作用處的值。3.1.5等效節(jié)點載荷單元等效節(jié)點力、和的表達式分別為

以上三式不但適用三角形單元，同樣適用其它單元，具有普遍意義。下面根據(jù)以上三式，分別進行討論。3.1.5等效節(jié)點載荷可以看出，當(dāng)外力作用在節(jié)點上時，如作用在i節(jié)點上，此時，，則計算所得到的等效節(jié)點力與原外力在數(shù)值上是相同的。3.1.5等效節(jié)點載荷（1）作用在單元任意點的集中力3.1.5等效節(jié)點載荷（2）作用在單元上的表面力現(xiàn)舉例說明，如圖所示的單元e，在ij邊上作用有表面力，假設(shè)ij邊的長度為l，單元厚度t為常數(shù)，其上任一點P距節(jié)點i的距離為s。有，，，則例3.3如圖所示，設(shè)31邊的長度為l，單元厚度為t，其上作用沿x方向的均布荷載，求其等效節(jié)點荷載。解：在31邊上，，

，則等效節(jié)點荷載為從計算結(jié)果來看，等效后的節(jié)點力是將原來作用在31邊上的均布面力平分到1、3節(jié)點上。3.1.5等效節(jié)點載荷（3）作用在單元上的體積力設(shè)單元所受體積力為，單元厚度t為常數(shù)，如所示，則3.1.5等效節(jié)點載荷當(dāng)單元體是均質(zhì)、等厚、比重為時，則、，有即當(dāng)單元有均勻自重時，在三節(jié)點三角形單元中每個節(jié)點沿y方向的等效節(jié)點力等于單元總重量的1/3。3.1.5等效節(jié)點載荷2.剛體靜力等效原則——力系簡化（1）作用在單元任意點的集中力若集中力作用在節(jié)點上，則可直接將集中力認為是等效節(jié)點力即可。若集中力作用于單元的某一邊界上，如圖所示，則先將分解為、，然后按線段的比例把分別移置到i，j兩點上。即3.1.5等效節(jié)點載荷（2）作用在單元上的表面力如圖所示，已知在ij邊受有均布面力，單元厚度t為常數(shù)，則移置到i、j節(jié)點上的等效節(jié)點力為如圖所示，當(dāng)某一邊上有三角形分布的面力時，則等效節(jié)點力為

3.1.5等效節(jié)點載荷（3）作用在單元上的體積力如圖所示，如果三角形單元ijm的重心c上受有體積力，單元厚度t為常數(shù)，則按剛體靜力等效原理可把體積力直接移置到i，j，m三個節(jié)點上，即根據(jù)圣維南原理可知，按剛體靜力等效原則移置后的等效節(jié)點力所引起的應(yīng)力誤差只是局部的不會影響到整體，即物體內(nèi)部距實際外力作用處相當(dāng)遠的各點處應(yīng)力，與其外力的具體分布情況關(guān)系很小。在線性位移模式下，彈性體和剛體按各自的靜力等效原則移置外力的結(jié)果是一致的。3.1.5等效節(jié)點載荷

整體節(jié)點載荷列陣：由各節(jié)點所受等效節(jié)點力按節(jié)點號碼以從小到大的順序排列組成的列陣。等效節(jié)點力是由集中力、表面力和體積力共同移置構(gòu)成的，其中集中力包括直接作用在彈性體上的外力和邊界約束力，如支座反力。前文對單元體進行了分折，得到了單元剛度方程，但要解決問題，還必須進一步建立整個計算模型的整體剛度方程。完成這一步的關(guān)鍵，在于怎樣將單元的剛度矩陣和節(jié)點荷載列陣，分別“組裝”成整體剛度矩陣和整體節(jié)點荷載列陣。這里通過研究任意節(jié)點的平衡來建立整體剛度矩陣，該方法不但比較直觀、易懂，而且對怎樣編寫計算機程序是很有幫助的。為了研究整體剛度矩陣的組裝過程，先引入兩個概念。

整體節(jié)點位移列陣：由各節(jié)點位移按節(jié)點號碼以從小到大的順序排列組成的列陣。3.1.6整體剛度矩陣及其特點式中：，

3.1.6整體剛度矩陣及其特點不失一般性，僅考慮計算模型中有4個單元，如圖所示。四個單元的整體節(jié)點位移列陣為對每個單元都可以寫出相應(yīng)的單元剛度方程，即單元節(jié)點平衡方程。例如，對①號單元，有式中：——①號單元中第i（i=1,2,3）節(jié)點所受力。為了便于組裝整體剛度矩陣，將上式以整體節(jié)點位移表示，即

3.1.6整體剛度矩陣及其特點——①號單元的擴大剛度矩陣或稱為單元貢獻矩陣。同理，對于②單元，有式中：——②號單元中第i（i=1,3,4）節(jié)點所受力；——②號單元的擴大剛度矩陣。3.1.6整體剛度矩陣及其特點對于④單元，有式中：——④號單元中第i（i=3,4,5）節(jié)點所受力；——④號單元的擴大剛度矩陣。3.1.6整體剛度矩陣及其特點對于任意一個節(jié)點，可能承受兩種力的作用，一種是其它單元給予該節(jié)點的反作用力；另一種是作用在節(jié)點上的等效節(jié)點力。對整體而言，前者屬于內(nèi)力，后者屬于外力，每個節(jié)點在兩種力的作用下處于平衡。將各單元剛度方程左邊相加，即將各節(jié)點所受力相加，由于對于整體而言，單元給予節(jié)點的反作用力屬于內(nèi)力，在相加過程中相互抵消，所以各節(jié)點所受力相加的結(jié)果只有外力，即等效節(jié)點力，從而得到整體節(jié)點荷載列陣，如下3.1.6整體剛度矩陣及其特點將各單元剛度方程右邊相加，從而得到整體剛度矩陣，如下3.1.6整體剛度矩陣及其特點通過以上分析得，整體節(jié)點載荷與整體節(jié)點位移之間的關(guān)系式，即結(jié)構(gòu)整體有限元方程，如下式中：——整體剛度矩陣。3.1.6整體剛度矩陣及其特點整體剛度矩陣組裝的基本步驟：1）將單元剛度矩陣中的每個子塊放到在整體剛度矩陣中的對應(yīng)位置上，得到單元的擴大剛度矩陣。注意對于單元剛度矩陣是按照局部編碼排列的，即對應(yīng)單元剛度矩陣中的i、j、m；對于整體剛度矩陣是按照整體編碼排列的，即按節(jié)點號碼以從小到大的順序排列。在組裝過程中，必須知道單元節(jié)點的局部編碼與該節(jié)點在整體結(jié)構(gòu)中的整體編碼之間的關(guān)系，才能得到單元剛度矩陣中的每個子塊在整體剛度矩陣中的位置。將單元剛度矩陣中的每個子塊按總體編碼順序重新排列后，可以得到單元的擴大矩陣。例如在圖中，單元②的局部編碼為i、j、m，對應(yīng)整體編碼為1、3、4，然后將單元②剛度矩陣中的每個子塊按總體編碼順序重新排列后，可以得到單元的擴大矩陣。注意有些書籍中將局部編碼表示為1、2、3或1，2，3等；2）將全部單元的擴大矩陣相加得到整體剛度矩陣。3.1.6整體剛度矩陣及其特點通過以上組裝過程可以得到組裝整體剛度矩陣的一般規(guī)則：1）結(jié)構(gòu)中的等效節(jié)點力是相關(guān)單元結(jié)點力的疊加，整體剛度矩陣的子矩陣是相關(guān)單元的單元剛度矩陣子矩陣的集成；2）當(dāng)整體剛度矩陣中的子矩陣中r=s時，該節(jié)點（節(jié)點r或s）被哪幾個單元所共有，則就是這幾個單元的剛度矩陣中的子矩陣的相加。如應(yīng)該是單元①－④中對應(yīng)子矩陣的集成，即3.1.6整體剛度矩陣及其特點3.1.6整體剛度矩陣及其特點3）當(dāng)中時，若rs邊是組合體的內(nèi)邊，則就是共用該邊的兩相鄰單元剛度矩陣中的子矩陣的相加。如13邊為單元①和②的共用邊，則4）當(dāng)中r和s不同屬于任何單元時，則=0。如節(jié)點r＝1和s＝5不同屬于任何單元，此時=0。上述組裝基本步驟和規(guī)則具有普遍意義，對于不同類型、不同形式的單元，只是相應(yīng)節(jié)點的子矩陣的階數(shù)（節(jié)點自由度×節(jié)點自由度）可能不同，至于組裝整體剛度矩陣的規(guī)律仍是相同的。正是應(yīng)為有了這種組裝規(guī)律，使得有限元法能夠很方便地應(yīng)用電子計算機進行計算。例3.4如圖所示有限元模型，彈性模量為，厚度為，為簡化計算取，求整體剛度矩陣。3.1.6整體剛度矩陣及其特點單元編號①②③④整體編碼1、2、32、4、55、3、23、5、6局部編碼i、j、mi、j、mi、j、mi、j、m以整體編碼表示的單元剛度矩陣子塊解：該模型中共有6個節(jié)點，4個單元，各單元的信息如表3.1所示。表3.1各單元信息3.1.6整體剛度矩陣及其特點同例3.2類似的分析，得根據(jù)單元剛度矩陣的性質(zhì)，可知，由于幾何矩陣，所以有，整體剛度矩陣中的各子塊是對所有單元相應(yīng)的子塊求和得到的（實際只是對相關(guān)單元求和），其中各子塊矩陣均為2行×2列，整體剛度矩陣用子塊矩陣可以表示為3.1.6整體剛度矩陣及其特點上式中任意一子塊矩陣均為2行×2列，在計算過程中，無需將每個單元剛度矩陣進行擴大，只需判斷整體剛度矩陣子塊的下標(biāo)，然后利用組裝整體剛度矩陣的一般規(guī)則進行計算，如，由圖形可知，節(jié)點2由單元①、②和③所共有，則3.1.6整體剛度矩陣及其特點，由圖形可知，25邊為單元②和③的共用邊，則，由圖形可知，節(jié)點1、5不同屬于任何單元，則采用同樣的方法進行計算，得到整體剛度矩陣為3.1.6整體剛度矩陣及其特點3.1.6整體剛度矩陣及其特點4.是奇異矩陣，在排除剛體位移后，它是正定陣3.是稀疏矩陣，非零元素呈帶狀分布用有限元方法分析復(fù)雜工程問題時，節(jié)點的數(shù)目比較多，整體剛度矩陣的階數(shù)通常也是很高的。那么在進行計算時，如果存儲整體剛度矩陣的全部元素，將會浪費較大的資源、降低計算效率。如果根據(jù)整體剛度矩陣的特點進行編寫程序，可以大大節(jié)省資源、并提高計算效率。因此有必要了解和掌握整體剛度矩陣的特點，整體剛度矩陣具有以下幾個顯著的特點：1.是對稱矩陣2.中主對角元素總是正的3.1.6整體剛度矩陣及其特點3.1.6整體剛度矩陣及其特點1.是對稱矩陣由單元剛度矩陣的對稱性和整體剛度矩陣的組裝過程，可知整體剛度矩陣必為對稱矩陣。利用對稱性，在計算機編寫程序時，只存儲整體剛度矩陣上三角或下三角部分即可。2.中主對角元素總是正的例如，剛度矩陣中的元素表示節(jié)點2在x方向產(chǎn)生單位位移，而其它位移均為零時，在節(jié)點2的x方向上必須施加的力；表示節(jié)點2在y方向產(chǎn)生單位位移，而其它位移均為零時，在節(jié)點2的y方向上必須施加的力。很顯然在此情況下力的方向應(yīng)該與位移方向一致，故應(yīng)為正號。3.1.6整體剛度矩陣及其特點3.是稀疏矩陣，非零元素呈帶狀分布如果遵守一定的節(jié)點編號規(guī)則，就可使矩陣的非零元素都集中在主對角線附近呈帶狀。整體剛度矩陣中的子矩陣只有當(dāng)下標(biāo)s等于r或者s與r同屬于一個單元時才不為零，這就說明，在第r雙行中非零子矩陣的塊數(shù)，應(yīng)該等于節(jié)點r周圍直接相鄰的節(jié)點數(shù)目加1。可見，中元素一般都不是填滿的，是稀疏矩陣，且非零元素呈帶狀分布。以下圖所示的單元網(wǎng)格為例，其整體剛度矩陣中的非零子塊（每個子塊為2行2列）的分布情況如下圖所示，圖中陰影部分表示該子塊不為零，其它子塊部位均為零。3.1.6整體剛度矩陣及其特點3.1.6整體剛度矩陣及其特點顯然，帶狀剛度矩陣的帶寬取決于單元網(wǎng)格中相鄰節(jié)點號碼的最大差值D。把半個斜帶形區(qū)域中各行所具有的非零元素的最大個數(shù)叫做剛度矩陣的半帶寬（包括主對角元），用B表示，如下

B=2(D+1)通常的有限元程序，一般都利用剛度矩陣的對稱和稀疏帶狀的特點，在計算求解中，只存儲上半帶的元素，即所謂的半帶存儲。因此，在劃分完有限元網(wǎng)格進行節(jié)點編號時，要采用合理的編碼方式，使同一單元中相鄰兩節(jié)點的號碼差盡可能小，以便節(jié)省存儲空間、提高計算效率。3.1.6整體剛度矩陣及其特點對于同樣的有限元單元網(wǎng)格，按照圖（a）的結(jié)點編碼，最大的半帶寬為B=2×(10－4+1)=14；按照圖（b）的結(jié)點編碼，最大的半帶寬為B=2×(10－2+1)＝18；按照圖（c）的結(jié)點編碼，最大的半帶寬為B=2×(10－6+1)＝10。(a)(b)(c)3.1.6整體剛度矩陣及其特點4.是奇異矩陣，在排除剛體位移后，它是正定陣無約束的彈性體（或結(jié)構(gòu)物）的整體剛度矩陣是奇異的，不存在逆矩陣，即關(guān)于位移的解不唯一。這是因為彈性體在外力的作用下處于平衡，外力的分量應(yīng)該滿足三個靜力平衡方程。這反映在整體剛度矩陣中就意味著存在三個線性相關(guān)的行或列，所以是個奇異陣，不存在逆矩陣。例如：設(shè)彈性體在外力的作用下處于平衡，這時相應(yīng)的解為，然后在給予彈性體以剛體位移而相應(yīng)的節(jié)點位移，這時，仍是問題的解，因為剛體位移不會破壞平衡。注：當(dāng)排除剛體位移后，整體剛度矩陣是正定矩陣。前文已經(jīng)提到在排除剛體位移后，整體剛度矩陣是正定的，方程才可求得唯一解。排除剛體位移可以通過引入邊界約束條件來實現(xiàn)，這里介紹兩種比較簡單的引入已知位移的方法。1．代入法2．乘大數(shù)法3.1.7邊界條件3.1.7邊界條件1．代入法該方法保持方程組仍為2n×2n系統(tǒng)，僅對整體剛度矩陣和整體載荷列陣進行修正。下面以一個只有四個方程的簡單例子加以說明，方程如下假定系統(tǒng)中節(jié)點位移、，則當(dāng)引入這些節(jié)點的已知位移之后，方程就變成3.1.7邊界條件若，則然后，用這組維數(shù)不變的方程來求解所有的節(jié)點位移。顯然，其解答仍為原方程的解答。在手算時，可直接將零位移約束所對應(yīng)的整體剛度矩陣中的行和列直接劃去，使得整體剛陣的維數(shù)變小，更便于手算。3.1.7邊界條件2．乘大數(shù)法將中與指定的節(jié)點位移相對應(yīng)的主對角元素乘上一個大數(shù)，同時將中的對應(yīng)元素換成指定的節(jié)點位移值、該大數(shù)與節(jié)點位移相對應(yīng)的主對角元素三者的乘積。若把此方法用于上面的例子，則方程就變成該方程組的第一個方程為解得，這種方法就是使中相應(yīng)行的修正項遠大于非修正項。3.1.7邊界條件1．代入法2．乘大數(shù)法在以上的兩種方法中，代入法接近人工解法，雖然該方法比較直觀，但該方法對剛度矩陣改變較多，程序效率不高。乘大數(shù)法對剛度矩陣改變較少，工作量較小，但相乘的“大數(shù)”若取得過大，求解時會發(fā)生“溢出”、若取得太小則會引起較大的誤差。根據(jù)前面的討論，現(xiàn)以三角形常應(yīng)變單元為例來說明應(yīng)用有限元法求解彈性力學(xué)平面問題的具體步驟。1）力學(xué)模型的確定。根據(jù)工程實際情況確定問題的力學(xué)模型，并按一定比例繪制結(jié)構(gòu)圖、注明尺寸、載荷和約束情況等；2）結(jié)構(gòu)進行離散化。將彈性體劃分為許多單元，并對節(jié)點進行編號，確定全部節(jié)點的坐標(biāo)值；對單元進行編號，并列出各單元三個節(jié)點的節(jié)點號；3）計算等效節(jié)點載荷，形成整體載荷列陣；4）計算單元剛度矩陣；5）組裝整體剛度矩陣；6）處理約束，引入邊界條件；7）求解線性方程組，得到節(jié)點位移和節(jié)點力；8）計算單元應(yīng)力；9）整理計算結(jié)果（后處理部分）。3.1.8算例例3.5如圖所示為一厚度t=1cm的均質(zhì)正方形薄板，邊長為2m，上下受均勻拉力=106N/m2，材料彈性模量為E，泊松比，不記自重，試用有限元法求其應(yīng)力分量。

3.1.8算例解：1）力學(xué)模型的確定。由于此結(jié)構(gòu)長、寬遠大于厚度，而載荷作用于板平面內(nèi)，且沿板厚均勻分布，故可按平面應(yīng)力問題處理，考慮到結(jié)構(gòu)和載荷的對稱性，可取結(jié)構(gòu)的1/4來研究。3.1.8算例表3.2節(jié)點坐標(biāo)值節(jié)點坐標(biāo)1234x0110y00112）結(jié)構(gòu)進行離散化。為了計算簡便并能說明整個計算過程，將該1/4結(jié)構(gòu)離散為兩個三角形單元。節(jié)點編號、單元劃分如圖所示，各節(jié)點的坐標(biāo)值如表3.2所示，各單元編碼如表3.3所示。3.1.8算例表3.3單元編碼單元號①②整體編碼1、2、33、4、1局部編碼i、j、mi、j、m3）計算等效節(jié)點載荷，形成整體載荷列陣。按照彈性體靜力等效原則，應(yīng)用公式（3.41），并參考例3.3可得單元②由均布載荷等效的節(jié)點力為3.1.8算例若按照剛體靜力等效原則，參考公式（3.45），同樣可得到以上結(jié)果。考慮到支座反力，整體載荷列陣為4）計算單元剛度矩陣。對于單元①有，局部編碼i、j、m對應(yīng)整體編碼1、2、3，利用公式（3.34）并參考例3.2，可得單元①的剛度矩陣為3.1.8算例對于單元②有，局部編碼i、j、m對應(yīng)整體編碼3、4、1，由單元剛度矩陣的性質(zhì)可知，單元②的剛度矩陣與單元①的剛度矩陣相同。

5）組裝整體剛度矩陣。在組裝整體剛度矩陣時，要注意單元局部編碼與整體編碼之間的關(guān)系，按照組裝整體剛度矩陣的一般規(guī)則和基本步驟，可得整體剛度矩陣為3.1.8算例

6）處理約束，引入邊界條件。由以上分析可知，整體有限元方程為3.1.8算例根據(jù)邊界約束條件：，采用代入法引入邊界條件，劃去整體剛度矩陣中1、2、4、7的行和列，得3.1.8算例

7）求解線性方程組，得到節(jié)點位移。求解上面方程組可得出相應(yīng)得節(jié)點位移為3.1.8算例所以，整體節(jié)點為列陣為代入整體有限元方程，便可得整體載荷列陣為若要求得每個單元所受的節(jié)點力，則可以從整體節(jié)點位移列陣中提取單元所對應(yīng)的節(jié)點位移，然后在利用，便可求得單元節(jié)點所受的力。如單元①所對應(yīng)的節(jié)點位移為3.1.8算例求得單元①所對應(yīng)的節(jié)點力為單元②所對應(yīng)的節(jié)點位移為求得單元②所對應(yīng)的節(jié)點力為注意，這里求得的單元節(jié)點位移、和單元所對應(yīng)的節(jié)點力、均對應(yīng)局部編碼。很顯然，每個單元均滿足力的平衡。8）計算單元應(yīng)力。先求出各單元的應(yīng)力矩陣、，然后再求得各單元的應(yīng)力分量，如下3.1.8算例9）整理計算結(jié)果（后處理部分）。對于三節(jié)點三角形單元，單元內(nèi)各點的應(yīng)力值相等，算出的應(yīng)力一般作為單元形心處的應(yīng)力。由于單元應(yīng)力為常數(shù)，整個結(jié)構(gòu)的應(yīng)力場呈階梯狀，在單元之間不連續(xù)。而工程上往往更加關(guān)心邊界和節(jié)點上的受力情況，因此，必須對所得到的應(yīng)力再次進行處理，得到更加合理的應(yīng)力場，并得到所需點上的應(yīng)力值。這里介紹兩種簡單的方法，一種方法稱為節(jié)點平均法，即把環(huán)繞某一節(jié)點的各單元的應(yīng)力加以平均作為該節(jié)點的應(yīng)力值。例如圖中節(jié)點3的應(yīng)力為3.1.8算例為了使通過這樣平均得來的應(yīng)力比較接近實際情況，要求環(huán)繞節(jié)點的各單元尺寸不應(yīng)相差太大。這種做法，對內(nèi)節(jié)點比較好，對邊界點則可能很差。因此，邊界節(jié)點處的應(yīng)力不宜直接由單元應(yīng)力平均來獲得，而要根據(jù)內(nèi)節(jié)點的應(yīng)力構(gòu)造插值函數(shù)推算出來。例如圖中邊界點1的應(yīng)力，可以先用節(jié)點平均法求得節(jié)點2、3、4處的應(yīng)力，在構(gòu)造相應(yīng)的插值函數(shù)推算邊界點1的應(yīng)力，如常用的拋物線插值公式如下3.1.8算例另一種方法稱為單元平均法，即把兩相鄰單元的應(yīng)力加以平均，用以表示公共邊界中點的應(yīng)力。為了由這樣平均所得到的應(yīng)力具有較好的精度，兩相鄰單元的面積不應(yīng)相差太大。如圖中單元②和③邊界的中點處的應(yīng)力為3.1.8算例在不同的有限元軟件中均具有各自的后處理方法，但無論后期怎樣處理，應(yīng)力場來源于應(yīng)力的計算，應(yīng)力的精度主要依賴于單元的尺寸、單元的類型（位移模式）。對于有限元這種數(shù)值計算方法，一般總是希望隨著網(wǎng)格的逐步細分所得到的解能夠收斂于問題的精確解。根據(jù)前面的分析，可知在有限元分析中，一旦確定了單元的形狀之后，位移模式的選擇將是非常關(guān)鍵的。由于載荷的移置、應(yīng)力矩陣和剛度矩陣的建立等等，都依賴于單元的位移模式，所以，如果所選擇的位移模式與真實的位移分布有很大的差別，那么就很難獲得良好的數(shù)值解。為了保證解答的收斂性，要求位移模式必須滿足以下三個條件，即1）位移模式必須包含單元的剛體位移。2）位移模式必須包含單元的常應(yīng)變。3）位移模式在單元內(nèi)要連續(xù)、且在相鄰單元之間的位移必須協(xié)調(diào)。3.1.9收斂準(zhǔn)則3.1.9收斂準(zhǔn)則

1）位移模式必須包含單元的剛體位移。也就是說，當(dāng)節(jié)點位移是由某個剛體位移所引起時，彈性體內(nèi)將不會產(chǎn)生應(yīng)變。所以，位移模式不但要具有描述單元本身形變的能力，而且還要具有描述由于其它單元形變而通過節(jié)點位移引起單元剛體位移的能力。2）位移模式必須包含單元的常應(yīng)變。每個單元的應(yīng)變一般都是包含著兩個部分：一部分是與該單元中各點的坐標(biāo)位置有關(guān)的應(yīng)變（即所謂各點的變應(yīng)變）；另一部分是與位置坐標(biāo)無關(guān)的應(yīng)變（即所謂的常應(yīng)變）。從物理意義上看，當(dāng)單元尺寸無限縮小時，每個單元中的應(yīng)變應(yīng)該趨于常量。因此，在位移模式中必須包含有這些常應(yīng)變，否則就不可能使數(shù)值解收斂于正確解。3.1.9收斂準(zhǔn)則3）位移模式在單元內(nèi)要連續(xù)、且在相鄰單元之間的位移必須協(xié)調(diào)。當(dāng)選擇多項式來構(gòu)成位移模式時，單元內(nèi)的連續(xù)性要求總是得到滿足的，單元間的位移協(xié)調(diào)性，就是要求單元之間既不會出現(xiàn)開裂也不會出現(xiàn)重疊的現(xiàn)象。通常，當(dāng)單元交界面上的位移取決于該交界面上節(jié)點的位移時，就可以保證位移的協(xié)調(diào)性。在有限單元法中，把能夠滿足條件1和2的單元，稱為完備單元，完備單元是收斂的必要條件。滿足條件3的單元，叫做協(xié)調(diào)單元或保續(xù)單元。前面討論過的三節(jié)點三角形單元，均能同時滿足上述三個條件，因此屬于完備的協(xié)調(diào)單元，完備的協(xié)調(diào)單元是收斂的充分條件。在某些梁、板及殼體分析中，要使單元滿足條件3比較困難，所以實踐中有時也出現(xiàn)一些只滿足條件1和2的單元，其收斂性往往也能夠令人滿意。特別是放松條件3的單元，即完備而不協(xié)調(diào)的單元，這時要保證其收斂必須通過分片試驗的測試，關(guān)于分片試驗的內(nèi)容請參見相關(guān)書籍。目前，完備而不協(xié)調(diào)的單元，已獲得了很多成功應(yīng)用。3.1.9收斂準(zhǔn)則在有限元分析中，將實際連續(xù)體分成許多單元體的組合后，根據(jù)線性或非線性位移的假定，人為地選擇一個位移場，通過這些措施所得到的模型比實際連續(xù)體的剛性要高。因而，近似模型的剛度是實際連續(xù)體剛度的上界。若選擇不協(xié)調(diào)單元，那么這種模型可能由于單元分離、疊加或單元之間形成鉸而降低剛性，這種影響就有可能使得不協(xié)調(diào)元比應(yīng)用協(xié)調(diào)元所得的結(jié)果要好。不過，應(yīng)用不協(xié)調(diào)單元事先不能肯定所得的剛度是真實剛度的上界。換句話說，不協(xié)調(diào)元不一定象協(xié)調(diào)單元那樣剛硬，可能比較柔軟，因此有可能會比協(xié)調(diào)單元收斂得快。3.1.9收斂準(zhǔn)則經(jīng)驗證明，根據(jù)巴斯卡（Pascal）三角形來選擇二維多項式的各項。在二維多項式中，如果包含有對稱軸一邊的某一項，那么就必須同時包含有另一邊的對稱項。選擇多項式位移模式時，還應(yīng)該要考慮的一個因素是，多項式中的項數(shù)必須等于或稍大于單元邊界上的外節(jié)點的自由度數(shù)。通常是取項數(shù)與單元的外節(jié)點的自由度數(shù)相等，取過多的項數(shù)是不恰當(dāng)?shù)?。從?jié)點平衡得到結(jié)構(gòu)整體有限元方程式的方法物理概念清晰，但這種方法在數(shù)學(xué)上不夠嚴(yán)謹(jǐn)。本節(jié)從勢能原理導(dǎo)出有限元方程，該方法是一種更加嚴(yán)謹(jǐn)、適用性更加廣泛的方法。由第2章可知，系統(tǒng)的總勢能為如果有集中力作用等式右邊還應(yīng)加一項集中力所作的功，如下3.1.10從能量原理推導(dǎo)剛度矩陣當(dāng)把連續(xù)體離散為有限個單元，并在各單元內(nèi)假設(shè)了位移模式，在單元內(nèi)應(yīng)力與應(yīng)變用節(jié)點位移（或單元中）來描述，、，在連續(xù)體上的整體積分等于在每個單元上分別積分再求和，如下由勢能原理可知，對于處于平衡狀態(tài)的彈性體，位移真實解應(yīng)總使總勢能取極小值，反之依然，則由，得3.1.10從能量原理推導(dǎo)剛度矩陣上式進一步可寫為如果有集中力作用，有上式進一步可寫為上式同式（3.55），即為整體結(jié)構(gòu)有限元方程。

3.1.10從能量原理推導(dǎo)剛度矩陣3.2矩形單元雖然三角形單元具有很好的“適應(yīng)性”，幾乎任何復(fù)雜邊界的彈性體總可以劃分為三角形，并且三角形單元計算公式簡單，但精度較低。三角形單元間雖然能夠保證位移連續(xù)，但應(yīng)力的精度較差，不能很好的反映彈性體內(nèi)應(yīng)力的準(zhǔn)確分布規(guī)律。為了提高計算精度，準(zhǔn)確反映彈性體內(nèi)的應(yīng)力狀態(tài)，可以采用一些較精密的單元類型。本節(jié)將介紹常用的矩形單元，它采用了比常應(yīng)變?nèi)切螁卧螖?shù)更高的位移模式，因而可以更好地反映彈性體中的位移狀態(tài)和應(yīng)力狀態(tài)。另外，對一些邊界比較規(guī)則且呈直線的平面結(jié)構(gòu)的分析，采用矩形單元較合適。這時單元總數(shù)可以減少，相應(yīng)的原始數(shù)據(jù)準(zhǔn)備工作和單元特征計算工作均可節(jié)省。3.2.1位移函數(shù)如圖所示的矩形單元，不失一般性，令矩形單元的長、寬分別為2a、2b。矩形單元有4個節(jié)點，共8個自由度，即共有8個節(jié)點位移，采用類似三角形單元的分析方法，同樣可以完成對矩形單元的力學(xué)特性分析。這里引入一個局部坐標(biāo)系、，這樣可以推出比較簡潔的結(jié)果。如圖所示，取矩形單元的形心o為局部坐標(biāo)系的原點，和軸分別與整體坐標(biāo)軸x和y平行，兩坐標(biāo)系存在有以下的坐標(biāo)變換關(guān)系式中：、——矩形形心處坐標(biāo)。矩形形心處坐標(biāo)以及矩形長、寬可由下式計算3.2.1位移函數(shù)在局部坐標(biāo)系中，節(jié)點i的坐標(biāo)是，其值分別為±1。如節(jié)點1在局部坐標(biāo)系下的坐標(biāo)為（－1，－1）。由于矩形有4個節(jié)點，共8個自由度，可以選擇有8個待定參數(shù)的位移模式，如下采用類似三角形的分析方法，將各節(jié)點的局部坐標(biāo)代入上式，便可以求得相應(yīng)的待定參數(shù)，然后將帶回式位移模式，可得到用節(jié)點位移表示的位移模式，如下3.2.1位移函數(shù)式中：——矩形單元的形函數(shù)，i＝1，2，3，4；——形函數(shù)矩陣；——單元節(jié)點位移列陣，，i＝1，2，3，4。3.2.1位移函數(shù)（i＝1，2，3，4）形函數(shù)的表達式為3.2.1位移函數(shù)引入符號，，i＝1，2，3，4，則上式可以統(tǒng)一寫為可以看出，矩形單元的形函數(shù)具有和三角形單元形函數(shù)同樣的性質(zhì)，即：形函數(shù)在各單元節(jié)點上的值，具有“本點是1、它點為零”的性質(zhì)；在單元內(nèi)任意點上，四個形函數(shù)之和等于1；單元任意一條邊上的形函數(shù)，僅與該邊的兩端節(jié)點坐標(biāo)有關(guān)。3.2.1位移函數(shù)有了單元的位移模式，就可以利用平面問題的幾何方程求出單元內(nèi)任意點的應(yīng)變，將位移代入幾何方程，得式中的應(yīng)變轉(zhuǎn)換矩陣的子塊（i＝1，2，3，4）為3.2.2應(yīng)變與應(yīng)力矩陣求得應(yīng)變之后，再將應(yīng)變代入物理方程，便可推導(dǎo)出以節(jié)點位移表示的應(yīng)力，如下式中，應(yīng)力矩陣為其子塊（i＝1，2，3，4）為3.2.2應(yīng)變與應(yīng)力矩陣由前面的討論可以發(fā)現(xiàn)，四邊形單元的位移模式，式（3.77）比常應(yīng)變?nèi)切螁卧捎玫木€性位移模式增添了項（即相當(dāng)于xy項），把這種位移模式稱為雙線性模式。在這種模式下，單元內(nèi)的應(yīng)變分量將不再是常量，這一點可以從的表達式中看出。另外四邊形單元的位移模式中的與三角形單元相同，它反映了剛體位移和常應(yīng)變，而且在單元的邊界上(=±1或

=±1)，位移是按線性變化的，顯然在兩個相鄰單元的公共邊界上，其位移是連續(xù)的。3.2.2應(yīng)變與應(yīng)力矩陣由單元的應(yīng)力矩陣表達式還可以看出，矩形單元中的應(yīng)力分量也都不是常量。正應(yīng)力、和剪應(yīng)力均沿、兩個方向線性變化，即沿x、y兩個方向線性變化。正因為如此，若在彈性體中采用相同數(shù)目的節(jié)點時，矩形單元的精度要比常應(yīng)變?nèi)切螁卧木雀?。但是，矩形單元也有一些明顯的缺點，矩形單元不能適應(yīng)斜交的邊界和曲線邊界，不便于對不同部位采用不同大小的單元，以便提高有限元分析計算的效率和精度。3.2.2應(yīng)變與應(yīng)力矩陣矩形單元剛度矩陣的推導(dǎo)過程與三節(jié)點三角形單元類似，同樣可以表示為式（3.30）的形式，即，由前文可知的推導(dǎo)過程與形函數(shù)的具體表達形式、節(jié)點個數(shù)均無關(guān)，該表達式具有普遍意義。若單元厚度t

為常量，則可以進一步表示為將單元剛度矩陣寫成子塊的形式，如下3.2.3單元剛度矩陣上式中每一個子塊矩陣均為2行×2列，單元剛度矩陣中的子塊矩陣的表達式為

（r、s＝1，2，3，4）將應(yīng)變轉(zhuǎn)換矩陣子塊和彈性矩陣，代入上式，得3.2.3單元剛度矩陣式中：（r、s＝1，2，3，4）3.2.3單元剛度矩陣?yán)?.6如圖3.22所示，該模型中有兩個四邊形單元，彈性模量為＝210GPa，厚度為＝0.025m，泊松比＝0.3，＝1kN，求單元所受應(yīng)力。3.2.4算例解：單元①所對應(yīng)的節(jié)點為2、4、3、1，單元②所對應(yīng)的節(jié)點為4、6、5、3。在有限元分析過程中，首先求解單元①和②的剛度矩陣；然后組裝整體剛度矩陣，組裝的過程同三角形單元，此處不再給出；最后引入邊界條件（，＝0）并結(jié)合受力情況，求得整體節(jié)點位移列陣＝10-5×{0，0，0，0，0.1162，-0.1674，-0.1149，-0.1628，0.1514，-0.4707，-0.1568，-0.4978}T。為了求解單元應(yīng)力，需要先求得每個單元所對應(yīng)的節(jié)點位移列陣，結(jié)合單元的節(jié)點編號，可從整體位移列陣中提取單元①和單元②的位移列陣，如下＝10-5×{0，0，-0.1149，-0.1628，0.1162，-0.1674，0，0}T＝10-5×{-0.1149，-0.1628，-0.1568，-0.4978，0.1514，-0.4707，0.1162，-0.1674}T3.2.4算例注意：整體節(jié)點位移列陣是按照節(jié)點編號由小到大排列的，而單元位移列陣是按照單元節(jié)點編號排列的，如單元①所對應(yīng)的節(jié)點為2、4、3、1，則單元①的位移列陣中的前兩個數(shù)則表示節(jié)點2的x和y方向位移。將單元①和單元②的位移列陣代入應(yīng)力矩陣，可求得單元①和單元②的應(yīng)力，如下3.2.4算例通過以上兩式便可以求得單元①和單元②內(nèi)任意點的應(yīng)力，以單元①為例，若取，，則表示單元①的13邊中點處的應(yīng)力；若取，，則表示24邊中點處的應(yīng)力。若計算單元形心處的應(yīng)力，則取，為通過分析結(jié)果可知，單元內(nèi)任意點的應(yīng)力是坐標(biāo)的函數(shù)，若在彈性體中采用相同數(shù)目的節(jié)點時，矩形單元的精度顯然要高于常應(yīng)變?nèi)切螁卧木取?.2.4算例3.3平面等參元前一節(jié)討論的矩形單元具有精度較高、形狀規(guī)則、便于計算自動化等優(yōu)點，對于結(jié)構(gòu)復(fù)雜的曲邊外形，只能通過減小單元尺寸，增加單元數(shù)量進行逐漸逼近。這樣，自由度的數(shù)目隨之增加，計算時間長，工作量大。另外，這些單元的位移模式是線性模式，是實際位移模式的最低級逼近形式，問題的求解精度受到限制。為了克服以上缺點，人們試圖找出這樣一種單元：一方面，單元能很好地適應(yīng)曲線邊界和曲面邊界，準(zhǔn)確地模擬結(jié)構(gòu)形狀；另一方面，這種單元要具有較高次的位移模式，能更好地反映結(jié)構(gòu)的復(fù)雜應(yīng)力分布情況，即使單元網(wǎng)格劃分比較稀疏，也可以得到比較好的計算精度。等參數(shù)單元（等參元）就具備了以上兩條優(yōu)點，等參元為任意四邊形單元，其網(wǎng)格劃分不受邊界形狀的限制，單元大小可以不相等，是一種精度高而且應(yīng)用廣泛的單元。3.3平面等參元然而直接對任意四邊形進行單元分析是困難的，這是由于它的幾何形狀不規(guī)則，沒有統(tǒng)一的形狀，對各個單元逐個按不同的公式計算，因其工作量大而難以進行。為了解決這一問題，人們采用坐標(biāo)變換的方法，如圖所示，把坐標(biāo)系ξη（局部坐標(biāo)系）內(nèi)形狀規(guī)則單元（母單元）變換為另一坐標(biāo)系xy（整體坐標(biāo)系）中的任意四邊形單元（子單元）。然后導(dǎo)出關(guān)于局部坐標(biāo)系形狀規(guī)整的單元（母單元）的高階位移模式的形函數(shù)，然后利用形函數(shù)進行坐標(biāo)變換，得到關(guān)于整體坐標(biāo)系的復(fù)雜形狀的單元（子單元）。a）母單元b）子單元3.3平面等參元從圖形變換的角度看，坐標(biāo)系ξη和坐標(biāo)系xy可以分別看成是母單元和子單元這兩個不同單元的坐標(biāo)系，它們都是直角坐標(biāo)系。而從另一角度看，坐標(biāo)系ξη和坐標(biāo)系xy又可以看成是同一單元（子單元）的兩種不同的坐標(biāo)系。坐標(biāo)系xy是子單元的直角坐標(biāo)系，而坐標(biāo)系ξη可看成是子單元的曲線坐標(biāo)系。可以看出坐標(biāo)系xy始終扮演同一角色，即子單元的直角坐標(biāo)；而坐標(biāo)系ξη則扮演兩種角色，它既是母單元的直角坐標(biāo)，又是子單元的曲線坐標(biāo)。在有限元分析中，兩者的作用是不同的。直角坐標(biāo)系xy在整個結(jié)構(gòu)的所有子單元中共同采用，所以稱為整體坐標(biāo)系。而曲線坐標(biāo)系ξη則只適用于單個獨立的子單元，所以稱為局部坐標(biāo)。整體坐標(biāo)在整體分析中采用，局部坐標(biāo)則在單元分析中采用。平面問題中常用的等參元有四節(jié)點四邊形單元、八節(jié)點曲邊四邊形單元和6~8可變節(jié)點曲邊四邊形單元等，本節(jié)以八節(jié)點曲邊四邊形等參元為例介紹平面等參元的分析過程。八節(jié)點曲邊四邊形等參元的母單元是二維二次單元。八個節(jié)點分別為正方形的四個角點和四個邊中點，母單元（邊長為2的正方形）采用直角坐標(biāo)系ξη。由于母單元為16個自由度，因此單元的位移模式為3.3.1等參元剛度矩陣式中：式中的16個常數(shù)用8個節(jié)點的位移（，）表示后，則上式為3.3.1等參元剛度矩陣采用坐標(biāo)變換可使母單元的八個節(jié)點坐標(biāo)與等參元的八個節(jié)點坐標(biāo)建立一一對應(yīng)的關(guān)系，整體坐標(biāo)和局部坐標(biāo)的變換式為式中：——母單元的形函數(shù)。可知，母單元上任意直邊，例如2－3邊，有ξ＝1，則有關(guān)的形函數(shù)N2、N3和N6均是η的二次函數(shù)，其余的形函數(shù)均為零。這樣變換就成為二次非線性變換，它可將母單元的直邊2－3映射成子單元的曲邊2－3。根據(jù)等參元的思想，利用上述的形函數(shù)，可得等參元的位移模式為3.3.1等參元剛度矩陣有了單元的位移模式，就可以利用平面問題的幾何方程求出單元內(nèi)任意點的應(yīng)變，將上式代入式（2.12），得等參元的應(yīng)變矩陣為式中：——i節(jié)點位移，（i=1，2，…，8）；3.3.1等參元剛度矩陣——單元幾何矩陣，其中為（i=1，2，…，8）在上式中，是ξ、η的函數(shù)，因此必須用坐標(biāo)變換式來轉(zhuǎn)換導(dǎo)數(shù)關(guān)系，根據(jù)復(fù)合函數(shù)的求導(dǎo)法則，有（i=1，2，…，8）3.3.1等參元剛度矩陣式中的J稱為雅可比（Jacobi）矩陣，為則3.3.1等參元剛度矩陣雅可比（Jacobi）矩陣的逆矩陣為將單元應(yīng)變代入平面問題的物理方程式，就得到平面八節(jié)點等參元的應(yīng)力列陣，為（i=1，2，…，8）3.3.1等參元剛度矩陣等參元剛度矩陣仍可用式（3.30）計算，即在上式中，應(yīng)變矩陣為局部坐標(biāo)系ξη下的顯式，因此在計算上式時應(yīng)將微面積dxdy變換成dξ和dη所包圍的面積dA，由于所以，有3.3.1等參元剛度矩陣則等參元剛度矩陣為應(yīng)該指出，上式是對ξ和η的重積分，盡管其積分區(qū)域十分簡單，但其被積函數(shù)卻比較復(fù)雜，一般該積分沒有顯式，需要采用數(shù)值積分法求解。至于等參元的等效節(jié)點力的計算一般也要采用數(shù)值積分的方法求解，這里不再討論，可參考相關(guān)書籍。由以上分析可知，在劃分單元時，只需確定單元節(jié)點的整體坐標(biāo)值，而不必畫出等參元的具體形狀，因為在計算中實際使用的只有單元八個節(jié)點在整體坐標(biāo)系下的位移值。等參元變換的條件為，因此在有限元網(wǎng)格劃分時，要特別注意這一點。3.3.1等參元剛度矩陣如圖所示，在圖b中3、4點退化為一個點，在該點，同理在圖c中2、3點退化為一個點，在該點，在圖d在點1、2、3處>0，在點4處<0，又因為在單元內(nèi)連續(xù)變換，所以單元內(nèi)肯定存在＝0，這是由于單元過分歪曲造成的。a）正常b）不正常c）不正常