1、用于振動分析的有限元方法指導老師：陳益 報告人：成志斌 韓宗彪 何瑜 寧鵬內容有限元介紹單個元素的運動方程整個系統的運動方程整個系統的邊界條件的加載及質量矩陣 MATLAB實例及總結單個元素的質量矩陣、剛度矩陣、力矢量及其轉化有限元法簡介 有限元法是一種可用于精確地(但近似)解決許多復雜的振動問題的數值方法。 對于基本的一維元素進行有限元分析，能得到質量矩陣與剛度矩陣和所需的力矢量，對于二維三維，元素矩陣會轉換成相關的更高維的空間。使用一致的和集中質量矩陣的有限元方程并結合邊界條件能為復雜系統提供解釋。 最后,使用MATLAB程序得到在軸向載荷下的指定節點位移,固有振動頻率和特征值分析。本章目

2、的 *認識用于解決不同類型振動問題的剛度和質量矩陣。 *將矩陣元素從局部坐標系變換到全球坐標系。 *裝配單元矩陣和應用邊界條件。 *對桿、梁元素進行靜態分析。 *對桿、梁元素進行動態分析來得到固有頻率和振型。 *在有限元振動分析使用一致的集中質量矩陣。 *使用MATLAB解決振動問題。有限元思想 1，實際結構被一些元素所取代,這些元素都是被假定為一個連續的結構部件即有限元，這些元素在特定點即節點上互相關聯。 2，如果解決方案的各方面都選擇得當,那么它可以收斂到精確的解決方案，因為組成總體結構的元素很小，在節點上的力的平衡和元素之間的位移都令人感到滿意,這樣整個結構(組合的元素)表現為單一實體。

3、 3，因為得到準確解很難，所以得到一個方便且逼近的近似解很有價值。元素的運動方程龍門刨銑床有限元模型三角板元素梁元素元素的運動方程 位移函數 形狀函數 各點對應位移 未知節點位移數 n 動能 應變能質量矩陣剛度矩陣 位移函數 形狀函數 各點對應位移 未知節點位移數 n 動能 應變能主要內容： 一，單元的質量、剛度矩陣、等效節點力 二，單元矩陣的坐標變換 三，整個系統的運動方程一 單元的質量、剛度矩陣，等效節點力矢量一 桿單元一個桿單元是從桿上劃分出的一個小段，如下圖所示。由于單元很小，、A均視為常量。現在就以這最簡單的桿單元，推導出它的質量、剛度矩陣，等效節點力。圖12.1圖12.1 (1)

4、求桿單元上任意點的位移u(x,t)本來，桿單元上任意點的位移u(x,t)與節點的位移u1(t)、u2(t)之間的關系是未知的，但是，只要單元劃分的足夠小，那么其間的關系就無關大局。所以可以假定它們之間有簡單的線性關系，即根據節點位移對單元內任意點位移進行插值：)()(),(2211tututxu（1）式中，1、 2稱為線性系數，與單元里點的位置有關，是x的函數。此函數與單元的形狀有關，又叫形狀函數。形狀函數和插值函數一樣是任意的，但必須邊界條件：)(), 0(1tutu)(),(2tutlu只有滿足此條件單元才能協調一致運動，而不致破壞系統的完整性，因此這兩個條件實際上就是變形協調條件。將式（

5、1）帶入（2）中，就可以得到形狀函數1(x)、2(x)所滿足的邊界條件： （2） 1 , 00 0 , 102121ll （3） 以上邊界條件確定了 xx21、 lxxlxx21 ,1由于這兩個函數的任意性我們可以用簡單的線性函數來近似，因此有：（4）代回（1）式中有： tulxtulxtxu21)1 (),(（5）為此，我們已經找到了用節點位移表示單元內任意一點位移的表達式。 （2）計算此單元的動能和勢能桿單元的動能可表示成：（6） 上式中，是材料的密度，A是桿單元的橫截面積。用矩陣形式表示（6）式為：（7）其中，（8）所以，質量矩陣可以認為是：（9） 桿單元的勢能可以寫成：（10）式中，E

6、是彈性模量，（10）表示成矩陣形式為：這里，所以剛度矩陣k可以表示成：（11）（12）(3)計算等效節點力設單元上x處作用有分布力f ( x , t)，現在要把它等效成節點力 tftf21,遵循等效原則，即原載荷和等效之后的節點載荷在虛位移上所做的虛功相等。其實，就是對應于廣義坐標 tutu21,的廣義力，為此，計算txf,所做的虛功：把上式寫成矩陣形式：（13）所以等效節點力可以寫成： tftuktum （14）二 梁單元 ,31tftf如下圖所示，一個梁單元也是有兩個節點，但是有四個自由度，每個節點處，有兩種位移形式，一個是線位移，即撓度，一種是角位移。圖中，tftf42, 是力，是力矩。

7、txf,是分布載荷 ,31tt tt42, 是對應的線位移，是對應的轉角。tx,是梁單元上任意位移 x處的撓度。圖12.2在靜載彎曲條件下，梁單元上任意點出的撓度是x的三次方程，可寫成：此方程必須滿足下面的邊界條件：由此可以求解處a (t)、b (t)、c (t)、d (t)，進而撓度方程為：（16）（17）（18）tan撓度的斜率上式可以寫成形狀函數的表示：其中，形函數分別為：梁單元的動能、勢能、虛功表達式分別為：（19）式中I是橫截面的慣性矩上式中：（20）（21）（22）通過上式，可以得到梁單元的質量、剛度矩陣，等效節點力：二 單元矩陣的坐標變換局部坐標系：以各個單元本身的軸線為基準所設

8、立的坐標系。便于計算節點位移。缺點：如果整個系統里各個單元取向各異，各個節點位移方向不一致，如下圖。如何使匯交于一個節點的各個桿件的節點位移真正相等?解決方法：進行坐標變換 如右圖的系統，有四個桿件， u 1(t) 、 u 2(t) 為局部坐標系的節點位移，U i 為全局坐標系下的位移圖12.3如下圖，節點位移在局部、全局坐標系中的關系：坐標變換矩陣（23）其中，因為單元的動能、勢能與坐標系無關:(24) 得到在全局坐標下的單元質量、剛度矩陣為：類似地，根據單元在兩個坐標系下的力所做的虛功相等：得到在全局坐標系下的等效節點力：三 全系統運動方程 經過坐標變換，各個單元的節點位移方向被統一起來，

9、但是不同的節點有不同的節點位移，為了便于綜合出全系統的運動方程，首先要建立全系統的節點位移向量。每個單元的節點位移向量與全系統的節點位移向量之間的關系：長方形矩陣由1和0組成單元節點位移向量單元節點位移向量 例如，圖的12.5中的1單元，方程，變為：把每個單元的動能相加，就得到了整個系統的動能：（把整個系統的動能表示成關于節點速度矢量的形式） 這樣就得到了整個系統的質量矩陣:類似地，考慮整個系統的勢能，便可以得到整個系統的剛度矩陣：整個系統的廣義力向量：最后得到整個系統的運動方程：12.6添加邊界條件 前文中，節點沒有固定，結構在節點力的作用下會發生剛體位移。 也就是說，矩陣K是奇異矩陣。通常

10、情況下，我們希望結構 的位移為零。 因此，我們需要添加邊界條件對矩陣M、K和向量F進行約束。 N：結構中自由節點位移的數目例1：桿件分析 如圖：均質；長0.5m；斷面截面積5e-4m2；楊氏模量200GPa；密度7850Kg/m3；左端固定。 a.節點2處施加1000N靜態軸向外力u2，求應力 b.求系統固有頻率 解： a： 平衡方程： A=5e-4，E=2e11，l = 0.5，f2=1000，代入方程得： u1：位移，f1：節點1處應力，添加邊界條件：u1=0，解得：u2=5e-10m 由應力 與應變 的關系： 表示長度的變化， 表示應變； b： 由剛度矩陣和質量矩陣，得特征值方程： 21

11、luu ll 式中w位固有頻率，U1、U2分別是節點1、2的振幅，添加邊界條件：U1=0；解得： 例2：梁的自然頻率解：梁被理想化為單一單元，局部和整體的節點位移相同，如 圖所示： 梁的剛度矩陣： 質量矩陣： 節點位移向量： 與端點相關的邊界條件：W1=0，W3=0；解得： 求解特征值： 乘以l/2EI得： 令系數矩陣的行列式等于0得： 方程的根即梁的自然頻率： 結果可以和精確解比較：12.7一致、集中質量矩陣 12.3節中推出的質量矩陣是一致質量矩陣，因為用于推導剛度矩陣的位移模型也用于推導質量矩陣。 一些動態問題可以用形式簡單的質量矩陣求解。最簡單的質量矩陣集中質量矩陣，可以通過將質點指定

12、到節點上。 集中質量針對平移和旋轉的元素, 假設在平均位置兩側的特定位移表現得像個剛體而剩余的元素不參與運動。 因此這種假設不包括元素位移之間存在的動態耦合,因此產生的元素質量矩陣是純粹的對角矩陣。 1、桿的集中質量矩陣： 2、梁的集中質量矩陣： 旋轉自由度的慣性影響被假定為0；若考慮慣性影響，有轉動慣量： 集中質量矩陣變為： 對于一般的動態問題，兩者誰能得到更精確的解？ 兩個質量矩陣很相似， 他們不考慮各種位移自由度的元素之間的動態耦合。 他們的形狀函數也近似，都是用靜態的位移模型推導而來。 然而,由于集中質量矩陣對角，在計算時他使用更少的存儲空間。 下面的例子說明了在一個簡單的振動問題中，

13、集中和一致質量矩陣的應用。 集中質量矩陣與一致質量矩陣： 例：桿的一致和集中質量矩陣 用一致和集中質量矩陣求如圖所示兩端固定桿的固有頻率，用兩個桿單元建模。解：單元的剛度和質量矩陣分別是： 質量矩陣的下標c和l分別表示一致和 集中質量矩陣。 由于該桿由兩個單元建模，組合的剛度和質量矩陣如下： 方框中的部分分別與單元1和2相關。 添加邊界條件U1=U3=0后，特征值問題為： 特征值w由以下方程課解： 代入已知條件得：用一致質量矩陣用集中質量矩陣 事實上，方程的精確解析解是： 解得：12.8 MATLAB應用舉例 例12.5 階梯軸的有限元分析 圖12.11中的階梯軸滿足一下條件：A1=1610-

14、4 m2 ， A2=910-4 m2 ， A3=410-4 m2 ，Ei=201010Pa，i=1,2,3，pi=7.8103Kg/m3， i=1,2,3，l1=1m，l2=0.5m，l3=0.25m。編寫一個MATLAB程序解決以下問題。 a，在在載荷p3=1000N下u1，u2，u3的位移 b，階梯軸的固有頻率和模態12.8 MATLAB應用舉例 解決方案：階梯軸的剛度矩陣和質量矩陣如下所示： 在載荷p3的作用下系統的平衡方程如下所示：解得：b.由下面方程可求得特征值K已在（E.4）中得到，M可由下式得到：由MATLAB求解（E.3）和（E.5）的程序如下所示例12.6 階梯梁的特征值分析的程序如圖12.12所示的兩端固定的階梯梁，編寫一個MATLAB程序，命名為Program17.m，對它進行特征值分析。Program17.m程序的輸入數據如下：xl(i)=階梯梁的長度 ixi(i)=階梯梁的轉動慣量 ia(i)=階梯梁的橫截面積 ib(i，j)=對應本地的自由度 j 和階梯梁的自由度 ie