【文章內(nèi)容簡介】 27） 比較方程（ 217）與（ 218），可知（ 219）中的位移 x 是（ 220）中復(fù)數(shù) x 的虛部，因此（ 225）的虛部就是方程（ 212）的特解，即有 )sin( ?? ?? tBx （ 228） 其中 B 為振幅， ? 為相位差。由式（ 226）、 223）及（ 224）得出穩(wěn)態(tài)強迫振動有如下的基本特點： 遼寧工程技術(shù)大學(xué)畢業(yè)設(shè)計（論文） 13 1)線性系統(tǒng)對簡諧激勵的穩(wěn)態(tài)響應(yīng)是頻率等同于激勵頻率而相位滯后與激振力的簡諧振動； 2)穩(wěn)態(tài)響應(yīng)的振幅及相位差只取決于系統(tǒng)本身的物理性質(zhì)和激振的頻率及力幅，而與系統(tǒng)本身進入運動的方式無關(guān)。 無阻尼系統(tǒng)對簡諧激勵的穩(wěn)態(tài)響應(yīng)可以從式（ 226）得出。 當(dāng) n??? 時，得到 1?? ， 0?? ，這時 tkPx ?? s in1 1 20 ?? （ 229） 當(dāng) n??? 時，得到 1?? ， ??? ，這時 )s in (1120 ??? ??? tkPx （ 223） 式（ 221）也可以寫成（ 222）的形式，這時相位差反映在振幅20 11??kP的符號中。上述結(jié)果也可以由直接設(shè) tBx ?sin? 并代入下列方程而得到： tPkxxm ?sin0???? （ 224） 為了具體討論影響穩(wěn)定響應(yīng)的振幅和相位差的各種因素，記 kPB 00? （ 225） 0B 實際是質(zhì)量塊在激振力幅靜作用下的最大位移。再引入無量綱的振幅放大因子? ，它定義為 2220 )2()1(1 ???? ???? BB （ 226） 由式（ 226）和（ 219）可以分別畫出以相對阻尼系數(shù) ? 為參數(shù)的曲線 —— ??? 曲線與 ??? 曲線，前者稱為幅頻響應(yīng)曲線 ，后者稱為相頻響應(yīng)曲線如圖所示 程序如下 for kesai=[,,] lamda=0::。 beta=1./(sqrt((1lamda.^2).^2+(2*kesai*lamda).^2))。 plot(lamda,beta) 王超：基于 MATLAB 的振動系統(tǒng)編程分析 14 hold on end axis([0 5 0 3])。 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 500 . 511 . 522 . 53頻率比振幅放大因子1 . 00 . 50 . 3 7 500 . 0 50 . 1 50 . 2 5 偏心質(zhì)量引起的強迫振動振幅放大因子 2222)2()1( ??? ? ???meMB (227) 程序如下： for kesai=[,,] lamda=0::。 beta=lamda./(sqrt((1lamda.^2).^2+(2*kesai*lamda).^2))。 plot(lamda,beta) hold on end axis([0 5 0 3])。 遼寧工程技術(shù)大學(xué)畢業(yè)設(shè)計（論文） 15 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 500 . 511 . 522 . 53頻率比MB/me0 . 50 . 1 00 . 1 50 . 2 50 . 3 7 50 . 5 01 . 0 支撐運動引起的強迫 振動振幅放大因子 2222)2()1( )2(1 ??? ??? ?? ??? aB (228) 程序如下： for kesai=[,,] lamda=0::。 beta=sqrt((1+(2*kesai*lamda).^2)./((1lamda.^2).^2+(2*kesai*lamda).^2))。 plot(lamda,beta) hold on end axis([0 5 0 3])。 王超：基于 MATLAB 的振動系統(tǒng)編程分析 16 0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 4 . 5 500 . 511 . 522 . 53頻率比振幅放大因子0 . 0 50 . 1 00 . 1 50 . 2 50 . 3 7 50 . 5 01 . 01 . 1 4 1 算例利用 MATLAB，繪制彈簧 質(zhì)量系統(tǒng)在簡諧力作用下的響應(yīng)曲線。已知數(shù)據(jù)如下： smxmxtNtFmNkkgm /,30c o s1 0 0)(,/2 0 0 0,5 00 ????? ?。 系統(tǒng)全解形式如下： tftfxtxtx nnnnn ???????? c o sc o s)(s i n)( 22 022 000 ?????? ? 式中， sr a dsr a dmkmFf n /30,/20,20510000 ?????? ?? 利用 MATLAB 繪制解曲線上式的程序如下： % F0=100。 wn=20。 m=5。 w=30。 遼寧工程技術(shù)大學(xué)畢業(yè)設(shè)計（論文） 17 x0=。 x0_dot=。 f_0=F0/m。 for i=1:101 t(i)=2*(i1)/100。 x(i)=x0_dot*sin(wn*t(i))/wn+(x0f_0/(wn^2w^2))*cos(wn*t(i)) +f_0/(wn^2w^2)*cos(w*t(i))。 end plot(t,x)。 xlabel(39。t39。)。 ylabel(39。x(t)39。)。 title(39。39。) 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 0 . 2 0 . 1 5 0 . 1 0 . 0 500 . 0 50 . 10 . 1 5tx(t)E x 3 . 1 1 本章小 結(jié) 基于 MATLAB 對單自由度自由振動繪制振動圖像，進行粘性阻尼，強迫振動振幅放大因子繪圖進行數(shù)據(jù)分析，使振動數(shù)據(jù)更加明顯。 王超：基于 MATLAB 的振動系統(tǒng)編程分析 18 3 基于 MATLAB 的多自由度系統(tǒng)編程分析 多自由度系統(tǒng) [16] 多自由度振動系統(tǒng)的數(shù)學(xué)模型： ? ?? ? ? ?? ? ? ?? ? ? ?M x C x K x f? ? ? （ 31） 其中 ? ?M 、 ??C 、 ??K 、 ??f 和 ??x 分別為質(zhì)量矩陣、阻尼矩陣、剛度矩陣、力向量和響應(yīng)向量。把這個時域矩陣方程變換到拉氏域（變數(shù)為 p ），并假定初始位移和初始速度為零，則得： ? ? ? ? ? ? ? ? ? ?2( ) ( ) ( )p M p C K X p F p? ? ? （ 32） 或 ? ?? ? ? ?( ) ( ) ( )Z p X p F p? （ 33） 式中 ? ?()Zp ：動剛度矩陣。 由式 32）或（ 33）可以得出傳遞函數(shù)矩陣 ? ?()Hp ： ? ? ? ?? ?( ) ( ) ( )X p H p F p? （ 34） 借助矩陣相關(guān)理論計算出來： ? ? ? ? ? ?1 ( ( ) )( ) ( ) ()a d j Z pH p Z p Zp??? （ 35） 式中 ? ?( ( ) )adj Z p ：為伴隨矩陣； ()Zp：為 ? ?()Zp 的行列式。 式 （ 35）的分 母，叫做系統(tǒng)的特征方程。類似單自由度系統(tǒng)，特征方程的根，即系統(tǒng)極點，決定系統(tǒng)的共振頻率。應(yīng)為 ??C 一般粘性阻尼矩陣，不一定滿足矩陣可對角化條件，為了把系統(tǒng)方程 （ 32） 轉(zhuǎn)化為一般特征值問題公式，需引入恒等式： ? ? ? ? ? ? ? ?( ) 0p M p M X?? （ 36） 將此 式與 （ 32） 合并： ? ? ? ? ? ? ? ?()p A B Y F ??? （ 37） 其中 遼寧工程技術(shù)大學(xué)畢業(yè)設(shè)計（論文） 19 ? ? ? ? ? ?? ? ? ?0 MAMC??? ????, ? ? ? ? ? ?? ? ? ?00MB K???? ????, ? ? ? ?? ?pXY X?????????, ? ? ? ?? ?0F F????? ??????? 令 ??F? =??0 , 則 （ 37） 的特征值滿足下列方程： ? ? ? ? 0p A B?? （ 38） 對于 N 自由度系統(tǒng)，此方程有 2N 個復(fù)共軛對出現(xiàn)的特征根： i i ii i ijj? ? ?? ? ???? ??? ?? ?? 其中 i? 阻尼因子； i? 為阻尼固有頻率。 第一階 固有頻率及主振型 [17,18] 在求解系統(tǒng)動力響應(yīng)時，系統(tǒng)較低的前幾階固有頻率及相應(yīng)的主振型占有重要的地位，為計算它們而采用下面的矩陣迭代法是比較簡單的。 將 i???? ?? 和i 帶入 公式 中，得 iii ??? ?A （ 39） 若將上式左端看作新列陣，上式表示：對于精確的主振型。新列陣 ( iA? ) 與原來的列陣 i? 的各個對應(yīng)元素之間都相差同一常倍數(shù)，這個常倍數(shù)即特征值 1? 。 記 1X 為初始迭代列陣，由展開定理， 1X 可以表示為 111 ?aX ? nnaa ?? ??? ?22 （ 310） 對上式左乘矩陣 A，由式（ 39）得知第一次迭代后所得的列陣為 nnnaaaAXX ?????? ????? ?22211112 = ???????? ??? nnnaaa ????????1212111?（ 311） 如果特征值 1? 不是特征方程的重根，那么上式中的1n1312 ?????? 、 ? 都小于 1，因此比起其他主振型 1? 在 2X 內(nèi)占的比重相對地比在 1X 中占的比重大，換句話說，用矩陣 A 迭代計算一次后，擴大了迭代列陣中第一階主振型的優(yōu)勢。經(jīng)第二次迭代后，得 王超：基于 MATLAB 的振動系統(tǒng)編程分析 20 ???????? ?????????????????????nnnaaaAXX ???????? 2122122112123 ? 同理第（ r1）次迭代后的結(jié)果為 ???????? ????????????????????? ???? nrnnrrrr aaaAXX ????????112112211111 ? （ 312） 可見隨著次數(shù)的增加，第一階主振型的優(yōu)勢越來越擴大，當(dāng)?shù)螖?shù)充分大時，由上式近似地得 1111 ?? aX rr ?? （ 313） 這時再迭代一次，得出 rrr XAXX 11 ???? （ 314） 由此看到迭代后的新列陣 1?rX 與原來列陣 rX 的各個對應(yīng)元素之間都僅相差一倍數(shù) 1? ，所以 rX 或 1?rX 就是對應(yīng)于 1? 的第一階主振型，而特征值 1? 可由下式算出 ? ?? ?lrlrXX 11 ??? )21( nl ，， ?? （ 315） 其中 ? ?表示lrX 列陣 rX 的第 l 個元素。為防止迭代過程中迭代列陣的元素變得過大或過小，每次迭代后需要使列陣歸一化，例如使它最后一個元素成為 1。下面是實用的矩陣迭代法的計算步驟： 1)選取初始迭代列陣 1X 使其最后一個元素為 1； 2)對 iX 作矩陣迭代，并使新 列陣 1Y 歸一化，即 11 AXY ? ， ? ? 1121 YYXn? （ 316） 3)重復(fù)步驟（ 2），第 r 次的迭代結(jié)果為 rr AXY ? ， ? ? rnrr YYX11??