用于振動(dòng)分析的有限元方法

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

2、的 *認(rèn)識(shí)用于解決不同類(lèi)型振動(dòng)問(wèn)題的剛度和質(zhì)量矩陣。 *將矩陣元素從局部坐標(biāo)系變換到全球坐標(biāo)系。 *裝配單元矩陣和應(yīng)用邊界條件。 *對(duì)桿、梁元素進(jìn)行靜態(tài)分析。 *對(duì)桿、梁元素進(jìn)行動(dòng)態(tài)分析來(lái)得到固有頻率和振型。 *在有限元振動(dòng)分析使用一致的集中質(zhì)量矩陣。 *使用MATLAB解決振動(dòng)問(wèn)題。有限元思想 1，實(shí)際結(jié)構(gòu)被一些元素所取代,這些元素都是被假定為一個(gè)連續(xù)的結(jié)構(gòu)部件即有限元，這些元素在特定點(diǎn)即節(jié)點(diǎn)上互相關(guān)聯(lián)。 2，如果解決方案的各方面都選擇得當(dāng),那么它可以收斂到精確的解決方案，因?yàn)榻M成總體結(jié)構(gòu)的元素很小，在節(jié)點(diǎn)上的力的平衡和元素之間的位移都令人感到滿(mǎn)意,這樣整個(gè)結(jié)構(gòu)(組合的元素)表現(xiàn)為單一實(shí)體。

3、 3，因?yàn)榈玫綔?zhǔn)確解很難，所以得到一個(gè)方便且逼近的近似解很有價(jià)值。元素的運(yùn)動(dòng)方程龍門(mén)刨銑床有限元模型三角板元素梁元素元素的運(yùn)動(dòng)方程 位移函數(shù) 形狀函數(shù) 各點(diǎn)對(duì)應(yīng)位移 未知節(jié)點(diǎn)位移數(shù) n 動(dòng)能 應(yīng)變能質(zhì)量矩陣剛度矩陣 位移函數(shù) 形狀函數(shù) 各點(diǎn)對(duì)應(yīng)位移 未知節(jié)點(diǎn)位移數(shù) n 動(dòng)能 應(yīng)變能主要內(nèi)容： 一，單元的質(zhì)量、剛度矩陣、等效節(jié)點(diǎn)力 二，單元矩陣的坐標(biāo)變換 三，整個(gè)系統(tǒng)的運(yùn)動(dòng)方程一 單元的質(zhì)量、剛度矩陣，等效節(jié)點(diǎn)力矢量一 桿單元一個(gè)桿單元是從桿上劃分出的一個(gè)小段，如下圖所示。由于單元很小，、A均視為常量。現(xiàn)在就以這最簡(jiǎn)單的桿單元，推導(dǎo)出它的質(zhì)量、剛度矩陣，等效節(jié)點(diǎn)力。圖12.1圖12.1 (1)

4、求桿單元上任意點(diǎn)的位移u(x,t)本來(lái)，桿單元上任意點(diǎn)的位移u(x,t)與節(jié)點(diǎn)的位移u1(t)、u2(t)之間的關(guān)系是未知的，但是，只要單元?jiǎng)澐值淖銐蛐?，那么其間的關(guān)系就無(wú)關(guān)大局。所以可以假定它們之間有簡(jiǎn)單的線(xiàn)性關(guān)系，即根據(jù)節(jié)點(diǎn)位移對(duì)單元內(nèi)任意點(diǎn)位移進(jìn)行插值：)()(),(2211tututxu（1）式中，1、 2稱(chēng)為線(xiàn)性系數(shù)，與單元里點(diǎn)的位置有關(guān)，是x的函數(shù)。此函數(shù)與單元的形狀有關(guān)，又叫形狀函數(shù)。形狀函數(shù)和插值函數(shù)一樣是任意的，但必須邊界條件：)(), 0(1tutu)(),(2tutlu只有滿(mǎn)足此條件單元才能協(xié)調(diào)一致運(yùn)動(dòng)，而不致破壞系統(tǒng)的完整性，因此這兩個(gè)條件實(shí)際上就是變形協(xié)調(diào)條件。將式（

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

6、是彈性模量，（10）表示成矩陣形式為：這里，所以剛度矩陣k可以表示成：（11）（12）(3)計(jì)算等效節(jié)點(diǎn)力設(shè)單元上x(chóng)處作用有分布力f ( x , t)，現(xiàn)在要把它等效成節(jié)點(diǎn)力 tftf21,遵循等效原則，即原載荷和等效之后的節(jié)點(diǎn)載荷在虛位移上所做的虛功相等。其實(shí)，就是對(duì)應(yīng)于廣義坐標(biāo) tutu21,的廣義力，為此，計(jì)算txf,所做的虛功：把上式寫(xiě)成矩陣形式：（13）所以等效節(jié)點(diǎn)力可以寫(xiě)成： tftuktum （14）二 梁?jiǎn)卧?,31tftf如下圖所示，一個(gè)梁?jiǎn)卧彩怯袃蓚€(gè)節(jié)點(diǎn)，但是有四個(gè)自由度，每個(gè)節(jié)點(diǎn)處，有兩種位移形式，一個(gè)是線(xiàn)位移，即撓度，一種是角位移。圖中，tftf42, 是力，是力矩。

7、txf,是分布載荷 ,31tt tt42, 是對(duì)應(yīng)的線(xiàn)位移，是對(duì)應(yīng)的轉(zhuǎn)角。tx,是梁?jiǎn)卧先我馕灰?x處的撓度。圖12.2在靜載彎曲條件下，梁?jiǎn)卧先我恻c(diǎn)出的撓度是x的三次方程，可寫(xiě)成：此方程必須滿(mǎn)足下面的邊界條件：由此可以求解處a (t)、b (t)、c (t)、d (t)，進(jìn)而撓度方程為：（16）（17）（18）tan撓度的斜率上式可以寫(xiě)成形狀函數(shù)的表示：其中，形函數(shù)分別為：梁?jiǎn)卧膭?dòng)能、勢(shì)能、虛功表達(dá)式分別為：（19）式中I是橫截面的慣性矩上式中：（20）（21）（22）通過(guò)上式，可以得到梁?jiǎn)卧馁|(zhì)量、剛度矩陣，等效節(jié)點(diǎn)力：二 單元矩陣的坐標(biāo)變換局部坐標(biāo)系：以各個(gè)單元本身的軸線(xiàn)為基準(zhǔn)所設(shè)

8、立的坐標(biāo)系。便于計(jì)算節(jié)點(diǎn)位移。缺點(diǎn)：如果整個(gè)系統(tǒng)里各個(gè)單元取向各異，各個(gè)節(jié)點(diǎn)位移方向不一致，如下圖。如何使匯交于一個(gè)節(jié)點(diǎn)的各個(gè)桿件的節(jié)點(diǎn)位移真正相等?解決方法：進(jìn)行坐標(biāo)變換 如右圖的系統(tǒng)，有四個(gè)桿件， u 1(t) 、 u 2(t) 為局部坐標(biāo)系的節(jié)點(diǎn)位移，U i 為全局坐標(biāo)系下的位移圖12.3如下圖，節(jié)點(diǎn)位移在局部、全局坐標(biāo)系中的關(guān)系：坐標(biāo)變換矩陣（23）其中，因?yàn)閱卧膭?dòng)能、勢(shì)能與坐標(biāo)系無(wú)關(guān):(24) 得到在全局坐標(biāo)下的單元質(zhì)量、剛度矩陣為：類(lèi)似地，根據(jù)單元在兩個(gè)坐標(biāo)系下的力所做的虛功相等：得到在全局坐標(biāo)系下的等效節(jié)點(diǎn)力：三 全系統(tǒng)運(yùn)動(dòng)方程 經(jīng)過(guò)坐標(biāo)變換，各個(gè)單元的節(jié)點(diǎn)位移方向被統(tǒng)一起來(lái)，

9、但是不同的節(jié)點(diǎn)有不同的節(jié)點(diǎn)位移，為了便于綜合出全系統(tǒng)的運(yùn)動(dòng)方程，首先要建立全系統(tǒng)的節(jié)點(diǎn)位移向量。每個(gè)單元的節(jié)點(diǎn)位移向量與全系統(tǒng)的節(jié)點(diǎn)位移向量之間的關(guān)系：長(zhǎng)方形矩陣由1和0組成單元節(jié)點(diǎn)位移向量單元節(jié)點(diǎn)位移向量 例如，圖的12.5中的1單元，方程，變?yōu)椋喊衙總€(gè)單元的動(dòng)能相加，就得到了整個(gè)系統(tǒng)的動(dòng)能：（把整個(gè)系統(tǒng)的動(dòng)能表示成關(guān)于節(jié)點(diǎn)速度矢量的形式） 這樣就得到了整個(gè)系統(tǒng)的質(zhì)量矩陣:類(lèi)似地，考慮整個(gè)系統(tǒng)的勢(shì)能，便可以得到整個(gè)系統(tǒng)的剛度矩陣：整個(gè)系統(tǒng)的廣義力向量：最后得到整個(gè)系統(tǒng)的運(yùn)動(dòng)方程：12.6添加邊界條件 前文中，節(jié)點(diǎn)沒(méi)有固定，結(jié)構(gòu)在節(jié)點(diǎn)力的作用下會(huì)發(fā)生剛體位移。 也就是說(shuō)，矩陣K是奇異矩陣。通常

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

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

12、到節(jié)點(diǎn)上。 集中質(zhì)量針對(duì)平移和旋轉(zhuǎn)的元素, 假設(shè)在平均位置兩側(cè)的特定位移表現(xiàn)得像個(gè)剛體而剩余的元素不參與運(yùn)動(dòng)。 因此這種假設(shè)不包括元素位移之間存在的動(dòng)態(tài)耦合,因此產(chǎn)生的元素質(zhì)量矩陣是純粹的對(duì)角矩陣。 1、桿的集中質(zhì)量矩陣： 2、梁的集中質(zhì)量矩陣： 旋轉(zhuǎn)自由度的慣性影響被假定為0；若考慮慣性影響，有轉(zhuǎn)動(dòng)慣量： 集中質(zhì)量矩陣變?yōu)椋?對(duì)于一般的動(dòng)態(tài)問(wèn)題，兩者誰(shuí)能得到更精確的解？ 兩個(gè)質(zhì)量矩陣很相似， 他們不考慮各種位移自由度的元素之間的動(dòng)態(tài)耦合。 他們的形狀函數(shù)也近似，都是用靜態(tài)的位移模型推導(dǎo)而來(lái)。 然而,由于集中質(zhì)量矩陣對(duì)角，在計(jì)算時(shí)他使用更少的存儲(chǔ)空間。 下面的例子說(shuō)明了在一個(gè)簡(jiǎn)單的振動(dòng)問(wèn)題中，

13、集中和一致質(zhì)量矩陣的應(yīng)用。 集中質(zhì)量矩陣與一致質(zhì)量矩陣： 例：桿的一致和集中質(zhì)量矩陣 用一致和集中質(zhì)量矩陣求如圖所示兩端固定桿的固有頻率，用兩個(gè)桿單元建模。解：?jiǎn)卧膭偠群唾|(zhì)量矩陣分別是： 質(zhì)量矩陣的下標(biāo)c和l分別表示一致和 集中質(zhì)量矩陣。 由于該桿由兩個(gè)單元建模，組合的剛度和質(zhì)量矩陣如下： 方框中的部分分別與單元1和2相關(guān)。 添加邊界條件U1=U3=0后，特征值問(wèn)題為： 特征值w由以下方程課解： 代入已知條件得：用一致質(zhì)量矩陣用集中質(zhì)量矩陣 事實(shí)上，方程的精確解析解是： 解得：12.8 MATLAB應(yīng)用舉例 例12.5 階梯軸的有限元分析 圖12.11中的階梯軸滿(mǎn)足一下條件：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。編寫(xiě)一個(gè)MATLAB程序解決以下問(wèn)題。 a，在在載荷p3=1000N下u1，u2，u3的位移 b，階梯軸的固有頻率和模態(tài)12.8 MATLAB應(yīng)用舉例 解決方案：階梯軸的剛度矩陣和質(zhì)量矩陣如下所示： 在載荷p3的作用下系統(tǒng)的平衡方程如下所示：解得：b.由下面方程可求得特征值K已在（E.4）中得到，M可由下式得到：由MATLAB求解（E.3）和（E.5）的程序如下所示例12.6 階梯梁的特征值分析的程序如圖12.12所示的兩端固定的階梯梁，編寫(xiě)一個(gè)MATLAB程序，命名為Program17.m，對(duì)它進(jìn)行特征值分析。Program17.m程序的輸入數(shù)據(jù)如下：xl(i)=階梯梁的長(zhǎng)度 ixi(i)=階梯梁的轉(zhuǎn)動(dòng)慣量 ia(i)=階梯梁的橫截面積 ib(i，j)=對(duì)應(yīng)本地的自由度 j 和階梯梁的自由度 ie