1、第第3章章 機械動力系統響應的數值計算機械動力系統響應的數值計算 本章主要內容：本章主要內容：3.1 歐拉法及其改進歐拉法及其改進3.2 線性加速度法線性加速度法3.3 紐馬克紐馬克- 法法3.4 威爾遜威爾遜法法3.4龍格龍格庫塔庫塔(RK)法法本章目的：本章目的：通過求解振動微分方程的通解，精確計算振動系統的響應通過求解振動微分方程的通解，精確計算振動系統的響應適適合求解單自由度系統。合求解單自由度系統。多自由度系統、非周期性激勵、非線性振動等：采用精確求解法多自由度系統、非周期性激勵、非線性振動等：采用精確求解法往往無法實現，只能數值計算方法。往往無法實現，只能數值計算方法。本章介紹的方

2、法，適用于單自由度、多自由度系統。對于多自由本章介紹的方法，適用于單自由度、多自由度系統。對于多自由度系統，各變量均用矩陣表示，度系統，各變量均用矩陣表示，M、C、K為方陣，為方陣，F、X均為列陣。均為列陣。 ，歐拉法算出的數值收斂于精確解。，歐拉法算出的數值收斂于精確解。 3.1 歐拉法及其改進歐拉法及其改進mxckxtfxttxtxttxttxtxttx/)()()()()()()( tttdttxtxttx)()()(dxtttxttxttxttt)()()(! 1)()( 0)(32)(!)()(! 3)()(! 2)()(! 1)()(kkktxkttxttxttxttxttx （3

3、.1）歐拉法：歐拉法：0t 誤差分析誤差分析: Taylor級數展開級數展開 特性特性:（3.2）（3.3）（3.4）構建位移表達式構建速度表達式加速度表達式動力學方程)(2)()()(2xtdxttRttt 討論討論（1）歐拉法是取）歐拉法是取Taylor級數展開式的前兩項的解法級數展開式的前兩項的解法（2）每前進一時間步）每前進一時間步 引起的誤差為引起的誤差為)(x ttt是是時間的加速度平均值時間的加速度平均值 其中其中 0 x 0 x （3）這意味著）這意味著R0時，歐拉法計算值過小，反之亦然時，歐拉法計算值過小，反之亦然波峰部分波峰部分，R0，即歐拉法計算值過小，即歐拉法計算值過小

4、與精確解相比，歐拉法數值解的振幅有逐漸增大的趨勢與精確解相比，歐拉法數值解的振幅有逐漸增大的趨勢改進的歐拉法改進的歐拉法（2）對速度表達式取具有）對速度表達式取具有更高精度的近似式更高精度的近似式 2/)()()()(2/)()()()(/)()()()(/)()()()(22ttxttxtxttxttxttxtxttxmtxctxktftxmtxctkxtftx tttxttxtxtxttx6)()2/(4)()()(（1）對）對Taylor級數級數取更高次項取更高次項，從而增加了計算量，從而增加了計算量( )x t提高了計算精度，但增加了提高了計算精度，但增加了tttxtxtxttx2)(

5、)()()(梯形法：梯形法： 辛普生辛普生(Simpson)公式：公式：速度表達式：Taylor級數取三項加速度表達式動力學方程三次加速度表達式位移表達式：Taylor級數取三項3.2 線性加速度法線性加速度法說明說明：式（3.9）、（）、（3.10）等號右邊含有）等號右邊含有 t+t 時刻量線性加速度法與歐拉法不同，它屬于隱式解法類型迭代計算。)(6)()(3)()()()(22ttxttxttx ttxttx 2)()()()(ttxtxttxttx ttxttxttxttx ttxttx)()(! 3)()(! 2)()()()(22 （3.10）)()()()(ttfttkxttxct

6、txm （3.11）（3.9）構建位移表達式：綜合梯形法構建位移表達式：綜合梯形法and辛普生辛普生(Simpson)公式公式構建速度表達式：梯形法構建速度表達式：梯形法大致相當于取到Taylor展開式的三次項。物理意義：假定從時刻tt+t時間的加速加速度直線變化度直線變化。聯立動力學方程：聯立動力學方程：線性加速度法直接解法一線性加速度法直接解法一FKXXCXM ktctmtxttx ttxktxttxcttfttx6/)()2/()(3)()()()(2)()()(22 （3.12）將（3.8）、 （3.9）代入（3.11），先消去對于多自由度系統)(6)()(3)()()()(2)()(

7、)()()(3)()()()(2)()(6)(2)(22212ttXttXttXttXttXttXtXttXttXtXttXttXKtXttXCttFKtCtMttX （3.13）和)(ttx)(ttx逆矩陣計算耗時，程序設計時，應置于循環之外2)()()()()(3)()()()()(6)()(6)()(12)()(3)()(2)(3)()()(6)(2)(22232212ttXtXttXttXtXttXttXttXtttXttFttXttXttXtCtXttXttXMKtCtMttX 線性加速度法直接解法二線性加速度法直接解法二)(ttx )(ttx222)/(6/3)()(3)(2)()

8、()(6)(6)(2)(tmtckttftxttxtxctxttxttxmttxt 將（3.8）、 （3.9）代入（3.11），先消去和FKXXCXM 對于多自由度系統多自由度系統（3.14）（3.16）逆矩陣計算耗時，程序設計時，應置于循環之外3.3 紐馬克紐馬克-法法 )()()()(2/)()()()(22tXttXttXttXttXttX 2/)()()()(ttXtXttXttX ttXttXttXttXttXttXttXtXtXttXtXttXttXKtXttXCttFKtCtMttx)()()()(! 2)()(! 1)()()()(2)()()(21()()()(2)()()(

9、2)(32212 ttxttxttxttx ttxttx)()(! 3)()(! 2)()()()(22 （3.10）上一節上一節“線性加速度法線性加速度法”中，構建的位移表達式中，構建的位移表達式構建位移表達式：構建位移表達式：構建速度表達式：構建速度表達式：梯形法（同前）梯形法（同前）（3.17）（3.18）紐馬克紐馬克-法法:紐馬克法是線性加紐馬克法是線性加速度法的別名。速度法的別名。 調節公式的特調節公式的特性參數，性參數，01/2。往往固定采用往往固定采用=1/6或或=1/4（3.19）3.4 威爾遜威爾遜法法 2)()()()()(6)()(3)()()()(22ttXtXttXt

10、tXttXttXttXttXttX 聯立動力學方程 )()()()(ttFttKXttXCttXM )(ttX )(ttX )()(1)(ttXtXttX ()()X ttX tt、6/ )()(3/ )()()()()(22ttXttXttXttXttX 2/)()()()(ttXtXttXttX 前面介紹的前面介紹的線性加速度法線性加速度法、紐馬克紐馬克-法法：均在時刻：均在時刻(t+t) 使用運使用運動方程動方程威爾遜威爾遜法法：應用于更后一點的時刻：應用于更后一點的時刻(t+ t )，1（3.20）先求，再用內插法求線性加速度法中線性加速度法中 t 替換成替換成 t而成而成再根據線性加

11、速度法構建式，求，即（3.21），即討論：威爾遜威爾遜法的物理意義：法的物理意義：加速度在時刻 t 到t+ t內為線性變化，首先計算t，t+ t 區間近似解，但僅取其中前半部分(到時刻t+ t為止)作為近似解，而舍去后半部分(時間t+ t以后)。這種巧妙的處理方法并非出于物理的原因，而主要是數學的理由。要理解這一點首先應了解數值計算的穩定性。振動仿真的失敗原因之一：振動仿真的失敗原因之一：往往是步長幅度過大。程序是正確的，輸入信息也對，但卻得到異常的結果。仔細研究可發現：計算剛開始時結果比較正常，但在計算過程中出現異?，F象，絕對值迅速增大。這種癥狀稱為不穩定。產生這種現象的原因很多，不一定只是

12、取值方面問題，但通常即使是良態方程，若 t過大，多半還是會出現不穩定現象。 t究竟取值：究竟取值：通常取小于周期的1/6。對多自由度系統，則 t應小于最短周期的1/6。在威爾遜法中，只要取大于1.37以上，不管取怎樣的值都是穩定的無條件穩定。因此，威爾遜法是實用價值很高的出色的解法，雖然增加了參數，式子稍復雜一些，但計算工作量與線性加速度法和紐馬克法差不多。威爾遜法兩種直接解法2122()() () ( )( )262() ( )( )( )3()(1)( )()/()( ) ( )()/2()( )( )()( )/3tttX ttMCKF ttC X tX ttK X ttX tX tX

13、ttX tX ttX ttX tt X tX ttX ttX ttX ttX t 2()()/3tX tt12223666() 2 ( )( )( )()()3( ) 2 ( )( )()266() ()( )( ) 2 ( )()11()(1) ( )()()( )CX ttKMMX tX tX ttttttCX tX tX tF tttX ttX ttX tX tX tttX ttX tX ttX ttX t 2(/2) ( )()()( )( ) () 2 ( )()/6tX tX ttX ttX ttX ttX tX tt 直接解法一直接解法二（常用）3.4龍格庫塔(RK)法（本科生略

14、） FxKxCxM 2200( ,)(0),(0)1,2,iiiiiiid xdxf t xdtdtxxxxin00( ,)(0),(0)iiiiiiiiidzf t x zdtdxzdtxxzx 對于多自由度振動系統（3-26）采用龍格庫塔(RK)法，既可以求解線性系統，也可以求解非線性系統n 維二階方程組降為一階方程方程組（3-27）（3-28）343342322312112114321143211),(2),2,2,2(2),2,2,2(),()22(6)22(6, 2 , 1hKzLhKzhLxhtfKKhzLKhzLhxhtfKKhzLKhzLhxhtfKzLzxtfKhLLLLhx

15、xKKKKhzzniiiiiiiiiiiiiiiiiiiii為步長四階龍格庫塔格式為（3-29）ODE（ordinary diffrential equation ）擴充內容：MATLAB中常微分方程解法介紹Matlab中求ODE數值解的函數Ode45四/五階龍格庫塔法，屬于單步法（只需要前一步的解就可以計算當前的解）。優點：不需要附加初始值，因此，計算過程中隨意改變步長也不會增加任何計算量。通常是很多問題首先試用的最好方法。缺點：不能求解剛性問題。Ode23四/五階龍格庫塔法，屬于單步法。在誤差允許范圍內較寬而且存在輕微剛性時，比Ode45效果好。Ode113可變階次的Adams-Bashforth-Moulton PECE算法。屬于多步法（需要前幾步的解來計算當前的解）。比Ode45更適合于誤差允許范圍要求嚴格的情況。缺點：不能求解剛性問題。Ode15s可變階數的NDFs算法。相對BDFs算法較好。是多步算法，剛性問題ode45不行時，可以試試這種算法。 Ode23t 自由內插方法的梯形算法，對剛性、又要求解沒有數值衰減時，可以用此法。 Ode23s改進的二價Rosenbrock算法。容許誤差較大時，ode23s比ode15好，所以