第3章 机械动力系统响应的数值计算-20121231

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以后)。这种巧妙的处理方法并非出于物理的原因，而主要是数学的理由。要理解这一点首先应了解数值计算的稳定性。振动仿真的失败原因之一：振动仿真的失败原因之一：往往是步长幅度过大。程序是正确的，输入信息也对，但却得到异常的结果。仔细研究可发现：计算刚开始时结果比较正常，但在计算过程中出现异常现象，绝对值迅速增大。这种症状称为不稳定。产生这种现象的原因很多，不一定只是

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好，所以