桿件結構有限單元法

2.1結構離散2.2單元的剛度矩陣2.3坐標變換2.4結構剛度方程2.5支座約束處理2.6剛度方程求解及內力計算例1.桁架結構計算示例例2.剛架結構計算示例第2章桿件結構分析的有限單元法2.1結構離散離散化要點：桿件的轉折點、匯交點、自由端、集中載荷作用點、支承點以及沿桿長截面突變處等均可設置成結點。結構中兩個結點間的每一個等截面直桿可以設置為一個單元。變截面桿件可分段處理成多個單元，仍按等截面桿單元計算，截面剛度取自各段中點處截面。對曲桿結構，可細分，用兩點之間的直線代替曲線。在有限元法計算中，載荷作用到結點上。當結構有非結點載荷作用時，應該按照靜力等效的原則將其變換為等效結點載荷。建立單元坐標系和結構整體坐標系。桿系結構離散化示意圖

2.1結構離散

2.1結構離散符號規(guī)定：線位移及相應力與坐標軸方向一致時為正；轉角位移和力矩，按右手法則定出的矢量方向若與坐標軸正向相一致時為正。基本未知量：鏈桿單元，平面每個節(jié)點兩個線位移，空間每個節(jié)點3個線位移梁單元，平面每個節(jié)點兩個線位移和一個轉角位移，空間每個節(jié)點3個線位移和三個轉角位移節(jié)點力與節(jié)點位移要一一對應平面剛架的梁單元e：平面桁架的鏈桿單元e：2.1結構離散6單元剛度方程反映單元的節(jié)點力和節(jié)點位移之間的關系;單元剛度矩陣可根據剛度系數的物理意義，由力和變形之間的關系確定。單元剛度矩陣也可以由位移函數和虛功原理或最小勢能原理計算。（具有一般性）2.2單元剛度矩陣返回7單元剛度矩陣（方法一）

對于一般平面剛架梁單元，單元剛度方程：單元的剛度方程：描述單元的結點位移δ(e)與結點力

F(e)之間的關系.結點位移結點力對于兩結點平面梁單元：單元的剛度矩陣K(e)為6*6的方陣對于平面桁架的鏈桿單元，K(e)為4*4的方陣。8平面剛架兩結點梁單元的剛度方程根據單元剛度系數的物理意義,由梁單元受力和變形及等截面直桿的剛度方程可以給出。 注意：所有力與位移的正負號均決定于坐標系方向。平面剛架兩節(jié)點梁單元：F(e)K(e)第i列元素的物理意義:

代表ui=1單獨作用于基本結構引起的六個桿端力大小。第krs元素的物理意義:

代表us=1單獨作用于基本結構引起的第r個桿端力Fr大小。δ(e)ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

ui=1平面梁單元的單元剛度矩陣單獨作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁單元的單元剛度矩陣vi=1單獨作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁單元的單元剛度矩陣單獨作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

uj=1平面梁單元的單元剛度矩陣單獨作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

vj=1平面梁單元的單元剛度矩陣單獨作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁單元的單元剛度矩陣單獨作用下15平面一般梁單元的單元剛度方程為:桿端力

向量F

(e)單元桿端位移

向量δ(e)單元剛度矩陣

K(e)16梁單元的單元剛度矩陣為:單元剛度矩陣常用子塊形式表示:其中每個都是3×3的方陣，子塊

Kij(e)

表示桿端j作用一單位位移時,桿i端引起的桿端力。返回17平面桁架鏈桿單元

平面桁架單元只有軸向變形,桿端力只有軸力；

矩陣表示：返回1）單元位移函數選用準則單元位移函數的項數，至少應等于單元的自由度數，它的階數至少包含常數項和一次項。位移函數中包含單元的剛體位移和常應變狀態(tài)。單元的位移函數應保證在單元內連續(xù)，以及相鄰單元之間的位移協(xié)調性。單元剛度矩陣（方法二）

根據結構分析類型選取單元類型和位移場，根據單元類型采用對應的單元位移模式，位移函數能反映真實結構的位移分布規(guī)律，保證計算精度及解的收斂性。選用的單元位移函數應當滿足下列要求:

形函數在其對應節(jié)點的值為1，在其他節(jié)點的值為0。2）軸向拉壓桿（鏈桿）單元的位移的函數由單元結點位移，代入位移函數中確定待定系數項：軸向位移形函數：軸向位移形函數：單元剛度矩陣（方法二）

梁單元撓度函數：3）梁單元平面彎曲的位移函數梁單元平面彎曲包括四個結點位移分量：由單元結點位移，代入位移函數中確定待定系數項：單元剛度矩陣（方法二）

稱為形函數矩陣。

梁單元撓度函數：梁單元位移模式：3）梁單元平面彎曲的位移函數單元剛度矩陣（方法二）

4）平面剛架梁單元的應力應變

平面剛架梁單元內任一點的軸向線應變由兩部分組成。則[B]—平面剛架梁單元的應變轉換矩陣。軸向應變：彎曲應變：平面剛架梁單元應變：

單元剛度矩陣（方法二）

5）平面剛架梁單元的剛度矩陣

梁單元的i，j結點發(fā)生虛位移為

單元內相應的虛應變應為：

由虛功原理有

由于結點虛位移的任意性，可得到單元的剛度方程：單元的剛度矩陣：單元剛度矩陣（方法二）

橫截面積

：橫截面對主慣性軸z的慣性矩：

通過積分同樣可以得到平面剛架梁單元的剛度矩陣：

平面剛架梁單元剛度矩陣25在整體分析時，要在結構整體坐標下進行桿端力（結點力）的疊加;通過坐標變換使所有單元的桿端力和桿端位移都變換到結構整體坐標下。在進行單元分析時，使用的是單元坐標系，各單元桿端力和桿端位移的排列順序和符號要參照單元坐標系。力和位移均為矢量，方向不同不能代數相加。2.3桿件單元的坐標變換返回設兩個坐標系之間的夾角為α，以整體到單元坐標系的轉向與x到y(tǒng)軸的轉向一致為正。26桿端力變換：桿端位移變換：平面鏈桿單元:平面桿件單元的坐標變換為正交矩陣:返回單元坐標變換矩陣：平面梁單元:整體坐標描述的單元剛度矩陣：

2.4結構剛度方程主要內容：231返回結構的剛度方程結構剛度矩陣的裝配結構的荷載向量結點位移向量Δ：所有結點的位移按順序排成一列；剛架每個結點都有3個位移，桁架每個結點都有2個位移。結點力向量P:作用在結點上的力按結點位移的順序排成一列;單元內部的荷載等效到節(jié)點上。返回結構剛度方程：反映結點力向量P與結點位移向量Δ之間的關系，即：KΔ=PK為結構的總剛度矩陣，且為對稱方陣。由個單元剛度矩陣裝配疊加而成。結點力向量P與結點位移Δ一一對應，也按結點位移的順序排列；以上所有的量都要用整體坐標表示。結構的剛度方程29具體做法：把整體坐標下的單元剛度矩陣根據結點編碼把各子塊送到總剛度矩陣K對應的位置中去。裝配過程：“子塊搬家，對號入座”。

整體坐標下單元剛度矩陣：如果i,j對應的結點編號為g,h，則單元剛度矩陣的各子塊在總剛中的位置分別為：結構剛度矩陣的裝配123ghn123ghn總剛度矩陣的集成子塊列子塊行每個單元的剛度矩陣都經過如上擴展和對號入座后，總剛度矩陣的各個子塊經過簡單的疊加即可得到最終的總剛度矩陣。“子塊搬家，對號入座”如：圖示平面剛架的總剛度矩陣的集成各單元結點編碼如圖。用矩陣記為：如果求出各單元整體坐標下的單元剛度矩陣：123456123456根據結點編號把各單元的子塊搬入總剛K中的對應位置。同一位置各子塊的對應元素相加；空位補0。如果修改各單元結點編號。則用矩陣：123456123456對角線下方劃紅色的元素需要存儲。單元的子塊搬入總剛度矩陣中的位置，完全取決于結構結點編碼。對同一結構，如果改變了結點的編碼，則總剛度矩陣完全不同。33總剛度矩陣的特點總剛度矩陣是一個對稱矩陣；處于主對角線對稱位置的兩個元素是相等的，即kij=kji

?？倓偠染仃囀且粋€稀疏的矩陣；大片的區(qū)域都是零元素，它的非零元素只分布在主對角線兩側的帶狀區(qū)域內。

最大半帶寬d=（c+1)*m;其中：c——各單元兩端結點編號差的最大值；

m——每個結點的自由度數；不相關結點對應的剛度子塊均為0?？倓偠染仃囀且粋€奇異矩陣；當沒有引進支座約束條件的情況下，總剛度矩陣不存在逆矩陣。返回判斷圖示結構總剛矩陣的最大半帶寬。d=(3+1)*3=1234d=(6+1)*3=21總剛矩陣中元素的排列與結點的順序直接相關。d=(7+1)*2=16返回最大相關結點差3最大相關結點差6最大相關結點差735結構荷載向量不考慮約束反力，只由外荷載引起的結點力排成的向量則稱為荷載向量。如圖1，結構的結點力向量:結構的荷載向量:荷載向量P的構成：直接結點荷載Pd等效結點荷載PE：單元內的非結點荷載（如分布荷載，溫荷載，慣性力等）等效移置到結點上得到的。

1234等效移置的方法:首先求出基本結構在非結點荷載作用下引起的固端力。最后將各固端力反向作用到結構的結點上去，即為該結點的等效結點荷載。1234等效結點荷載向量:荷載向量:2312荷載等效圖基本結構:兩端固定梁返回例4：計算圖示結構的荷載向量。

返回荷載等效圖

荷載等效圖

382.5支座約束處理在形成結構剛度方程時沒有考慮支座等約束條件;總剛度矩陣是一個奇異矩陣；剛度方程P=K沒有唯一解，方程中包含任意大小的剛體位移。必須引進約束條件，消除剛體位移，才能得到唯一解。支座條件引進的目的就是使：約束結點，位移已知0,約束反力未知。

自由結點，荷載已知,位移是未知量。結構剛度方程：K11(1)K12(1)0000K21(1)K22(1)+K22(2)0K24(4)0000K33(3)K34(2)0000K43(2)K44(2)+K44(3)+K44(4)0K46(4)0K52(4)00K55(5)K56(5)000K64(4)K65(5)K66(4)+K66(5)123456789101112131415161718P1P2P3P4P5P6P7P8P9P10P11P12P13P14P15P16P17P18對總剛度方程P=KΔ初等變換--行列交換K11(1)00K12(1)000K33(3)00K34(2)000K55(5)K52(4)0K56(5)K21(1)00K22(1)+K22(2)K24(4)00K43(2)00K44(2)+K44(3)+K44(4)K46(4)00K65(5)0K64(4)K66(4)+K66(5)123789131415456101112161718P1P2P3P7P8P9P13P14P15P4P5P6P10P11P12P16P17P18Px

未知支座反力P1已知節(jié)點

荷載Δ1支座已知0位移Δx求解的未知量Kx1KxxK11K1x把初等變換后的總剛度方程KΔ=

P可寫成將方程式展開得:已知節(jié)點力（節(jié)點荷載）未知節(jié)點力（支座反力）

未知節(jié)點位移求解的未知量已知節(jié)點位移（支座0位移)求解未知節(jié)點位移Δx計算未知約束反力Px對于剛性支座：支座反力對位移的計算沒有影響，但位移決定支座反力。有限元法是面向計算機的方法，邊界處理不要改變剛度方程的階數，消除剛體位移且不影響自由結點的位移。常用的處理方法有：1.主元素置1法：要修正剛度矩陣和荷載向量；作法：0位移對應的剛度矩陣主元素置為1，相應的副元素置0；荷載對應項置為0。2.主元素乘大數法：僅修正剛度矩陣；作法：僅把0位移對應的剛度矩陣主元素乘以一大數（如1030）。返回支座約束條件的引進主元素置1法若第i個自由度（位移）為0，應保證剛度方程解得Δi=0。

將總剛度矩陣K中第i行的主元素（第i行的主對角線元素）改為1，即令K（i，i）=1。將第i行、i列的所有副元素都改為零。即令K（i，j）=0，K（j，i）=0(i≠j)。將荷載向量中的第i元素置為零，即令Pi=0。

經過這三步改動后，便可實現Δi=0。

第3步第2步第1步第i行方程支座約束條件的引進123456789101112131415161718主元素置1法引進支座條件：

P1P2P3P4P5P6P7P8P9P10P11P12P13P14P15P16P17P18000P4P5P6000P10P11P12000P16P17P18根據0位移修改剛度矩陣；根據0位移修改荷載向量。(1)(2)(3)(4)(5)123456789101112131415161710P11P12P13P14P15P16P17P18如改變單元結點的編號；如圖引入邊界條件后的總剛方程為：主元素置1法引進支座條件：

返回形成了荷載向量P，集成了總剛度矩陣K并且引進支座條件后，便可由剛度方程K△=

P求解結點位移△。這就轉化為求解大型線性代數方程組問題。線性代數方程式組的解法：直接法和迭代法。直接解法：如高斯消去法，及其派生的LU、LDLT三角分解法。迭代法：塞德爾迭代法。2.6剛度方程求解及內力計算一、剛度方程求解返回47首先從整體的結點位移向量△中取出該單元的結點位移。設單元e左右兩端的結點編號分別為m和n，則該單元的整體坐標表示的結點位移向量為：然后將進行坐標變換，換成用單元坐標表示。最后代入單元剛度方程,便可求各單元的桿端力如果單元上還作用非結點荷載，則需要疊加由非結點荷載引起的固端內力，得到真正的桿端力：二、單元桿端內力計算返回例1.用有限單元法計算圖示桁架內力。EA=常量。

1.結構離散如圖所示2.單元坐標描述的各單元剛度矩陣結點坐標：1（0,1），2（1.732,1），3（0,0）；單元結點：1（1,2）,2（3,2）3.整體坐標描述的各單元剛度矩陣a=00a=300返回4.按結點編碼

裝配整體剛度矩陣單剛坐標變換：1231

2

300123212325.結構荷載向量6.根據邊界位移條件修正結構的剛度方程：5.結構荷載向量6.根據邊界位移條件修正結構的剛度方程：主元素置1法5.結構荷載向量6.根據邊界位移條件修正結構的剛度方程：主元素乘大數法7.求解結構的剛度方程組得位移：8.各單元的桿端位移向量：整體坐標下：轉換到單元坐標下：9.單元坐標下各單元的桿端力向量：單元內力：返回55例2.應用有限元法計算各桿端