计算动力学第1章课件

计算动力学（2）1§1.1振动的概念§1.2单自由度系统自由振动§1.3单自由度系统强迫振动§1.4两个自由度系统的振动§1.5非线性振动概述第1章绪论2§1.1振动的概念振动：就是物体在静平衡位置附近所作的往复运动。振动系统：在振动问题中所研究的对象。如机器或结构物等。激励：外界对振动系统的作用或引起机器运动的力。响应：机器或结构在激励作用下产生的动态行为。3力学基本模型振动系统的力学基本模型中包括三个基本“元件”：质量块、弹簧和阻尼器。质量块:是物体惯性大小的度量。弹簧:表示振动系统弹性的理想模型。阻尼器:任何振动在没有外界干扰(激励)时都会逐渐消失，因此，系统存在一种阻碍振动持续进行的阻力，这种阻力称为阻尼。5振动机理任何结构，之所以能产生振动，是因为它本身具有质量（惯性力）和弹簧（恢复力）。从能量关系看,质量可以储存动能,弹簧可以储存势能（变形能）。振动就是动能和势能不断地转换。6§1.2单自由度系统单自由度系统:可以用一个独立坐标来确定系统的位置及其运动规律的振动系统;单自由度线性系统的振动是最简单的振动系统;许多实际问题可以足够精确地简化为单自由度振动系统;单自由度振动系统的一些概念、特征和研究方法，是研究复杂振动系统的基础。7梁的横向振动质量为m的重物放在简支梁的中部，不计梁的质量。设梁长为l，材料的弹性模量为E，截面惯性矩为I。则利用材料力学的概念可得到：dst9m-c-k系统已知质量为m，弹簧的刚度系数为k，粘性阻尼系数为c。运动微分方程为：10m-c-k系统令阻尼比为则方程可写为令其解为代入方程得到此特征方程的两个根是11大阻尼情况给出初始条件：t＝0时则可确定系数B和D13大阻尼情况

这种情况对应的运动是一种衰减运动，但不是我们所关心的振动形式。设x0>0，v0>0，则运动图形大致如下。14临界阻尼情况（2）x＝1（临界阻尼情况）此时特征方程有重根利用初始条件确定常数为此时的阻尼系数称为临界阻尼系数，记为cc通解为15小阻尼情况（3）0<x<1（小阻尼情况）此时特征方程有一对共轭复根，通解为或写为利用初始条件确定出常数17小阻尼情况解中有两个因子，一个是衰减的指数函数,它将使振幅越来越小，直至振动最终消失;另一个是正弦函数，它表示系统以相同的周期通过平衡位置。18小阻尼情况因此系统呈现为一种衰减形式的等周期振动形式。19小阻尼情况（3）振动频率变小，周期变长此时系统振动的频率和周期为：因此：衰减振动的固有频率比无阻尼系统的固有频率小，振动周期变大，但影响不大，特别是当阻尼很小（x<<1）时，可以忽略阻尼对振动频率和周期的影响。21对数衰减率振幅衰减的快慢程度可用相邻振幅的比值来表示，称为衰减率或减幅率或减缩率；也可以用衰减率的自然对数来表示，称为对数衰减率。22对数衰减率利用前面给出的解可得到衰减率为对数衰减率为23对数衰减率则即对数衰减率则25§1.3单自由度系统强迫振动设激励为F(t)=F0sinwt，这里w为激振频率，利用牛顿定律并引入阻尼比x可得到26稳态响应性质2.幅频特性曲线对于稳态响应，定义动力放大系数R为响应的振幅X0与最大干扰力F0所引起的静位移的比值：以x为参数，画出R-r曲线即幅频特性曲线，表明了阻尼和激振频率对响应幅值的影响。29稳态响应性质Rr30§1.4两个自由度系统的振动单自由度系统振动问题，在我们所讨论的范围内是线性定常方程。而多自由度系统则是二阶多元联立微分方程组，各广义坐标间存在相互“耦合”现象。所谓耦合，就是变量之间互相联系。由于这种耦合，使微分方程的求解变得非常困难。因此，分析多自由度系统振动问题的重要内容之一就是如何将方程“解耦”，然后按单自由度的分析方法求解。两自由度是多自由度系统最简单的情况。31运动微分方程坐标原点仍取在静平衡位置写成矩阵形式32运动微分方程式中：33运动微分方程[M]称为系统的质量矩阵，[K]称为刚度矩阵，[C]称为阻尼矩阵，{x}为系统的位移列阵，{F(t)}为外激励列阵。对于其它形式的两自由度振动系统同样可得到相应的质量矩阵、刚度矩阵和阻尼矩阵。由于矩阵[M]、[K]、[C]的非对角线元素不为0，所以振动微分方程是互相耦合的非独立方程。34自由振动问题35自由振动问题特征方程特征根－纯虚根上述方程有非零解，要求系数矩阵的行列式为零36自由振动问题满足上述方程的特征向量37自由振动问题振型：第一阶振型第二阶振型方程的解38§1.5动力模型的建立1.拉格朗日方程2.有限元法39拉格朗日方程T系统动能V系统势能40拉格朗日方程例1：41拉格朗日方程例1(续）：42拉格朗日方程例1(续）：43拉格朗日方程例2：44拉格朗日方程例2(续)：45拉格朗日方程例2(续)：46拉格朗日方程例3：47拉格朗日方程例3(续)：48拉格朗日方程例3(续)：49拉格朗日方程例4：50拉格朗日方程例4(续)：静止状态：51拉格朗日方程例4(续)：运动状态：52拉格朗日方程例4(续)：53拉格朗日方程例4(续)：54两自由度系统例5：55两自由度系统例5(续)：56两自由度系统例5(续)：57三自由度系统例6：58三自由度系统例6(续)：动能：重力势能：当角度很小时：59三自由度系统例6(续)：弹簧伸长：弹簧势能：60三自由度系统例6(续)：总势能：61倒立单摆例7：62倒立单摆例7(续)：小车动能小球水平动能小球垂直动能小球势能外力势能63倒立单摆例7(续)：拉格朗日函数64倒立单摆例7(续)：65倒立单摆例7(续)：66倒立单摆例7(续)：67有限元拉压轴向压缩，对应的外力称为压力。轴向拉伸，对应的外力称为拉力。力学模型如图（截面积A，杆长L）68杆单元单位体积变形能69杆单元胡克定律单位体积变形能70杆单元单元总变形能71杆单元插值72杆单元形函数应变73杆单元节点力74杆单元刚度阵元素75杆单元刚度阵76梁单元弯曲力学模型如图（杆长L）zxM>0zxM<077梁单元单位体积变形能胡克定律78梁单元单位体积变形能单元总变形能79梁单元轴向位移应变80梁单元81梁单元横截面对y轴的惯性矩82梁单元插值83梁单元刚度阵84弹簧单元刚度阵85两自由度系统86刚度阵的建立87刚度阵的建立88§1.5非线性振动概述非线性特性

材料非线性振幅过大超出材料线弹性范围几何非线性位移或变形过大使结构几何形状显著变化非线性阻尼材料内摩擦阻尼、流体阻尼等都是非线性阻尼负刚度负阻尼有些情况下会存在负刚度和负阻尼非线性系统

当真实系统弹性元件的力与位移之间的关系超出线性范围，或阻尼元件的力与运动速度之间的关系不满足作线性关系时，系统的运动微分方程不能用线性微分方程描述，称系统为非线性系统。当真实系统作小运动时，可忽略系统的高阶微小量，近似地将系统看作线性系统。89非线性振动概述非线性振动的研究方法非线性振动研究的方法有：定性分析、定量分析和数值分析方法。非线性振动研究的内容非线性振动研究的基本内容之一就是建立对真实振动系统的计算方法，改进计算精度，探索某些特殊现象的规律。定性法研究已知解的领域内系统的一般稳定性特征，而不是运动的时间历程。通常采用几何方法描述系统的运动特征。定量法

通过一些渐近的解析方法研究系统运动的时间历程。数值法通过数值计算方法研究系统非线性振动的规律和现象。90非线性振动与线性振动的区别线性振动非线性振动

自由振动频率与初始条件无关自由振动频率与振幅有关

强迫振动频率与激励力频率相等强迫振动频率成分复杂，有时与激励频率不相等的频率成分突出稳定平衡位置附近的运动是稳定的稳定平衡位置附近具有多种稳定和不稳定运动强迫振动中每个激励频率有一个对应的振幅强迫振动中幅频与相频曲线发生弯曲，产生多值性

叠加原理成立叠加原理不成立91典型微分方程类型单摆方程库仑（Coulomb)摩擦振动方程92典型微分方程类型单摆方程库仑（Coulomb)摩擦振动方程93典型微分方程类型范德波（vanderPol）方程希尔(HiIl)方程94单摆

由牛顿第二定律：

非线性方程式中角频率：95单摆

线性化处理忽略3次以上的高次项得线性方程96单摆令代入方程得得特征方程：特征根：得通解为：式中为复常数。由于描述单摆振动的应为实函数，所以常数必须满足条件：将写成指数形式后得：该式是振幅为P，角频率为的简谐振动，其振动波形为正弦曲线。角频率只与摆线

l得长度有关，与摆锤质量无关，称为固有角频率。97单摆周期与摆角无关？看看实验结果：定性结论：1.周期随摆角增加而增加2.随摆角增加波形趋于矩形98单摆周期数学表达式对方程乘以后积分其中积分设t=0时，，周期为T，在时应有，故有：最后得：99§1.6常微分方程初值问题数值解法一阶常微分方程初值问题

在区间a≤x≤b上的数值解法。

(1)

100数值方法的基本思想101对常微分方程初值问题(1)式的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点处的函数值的近似值。相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定h为定数，称为定步长，这时节点可表示为101数值方法的基本思想102数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息 计算的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、102中的导数进行不同的离散化处理。对于初值问题数值方法的基本思想103数值微分、泰勒展开等离散化方法，对初值问题103的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算yi+1时只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格—库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法；其代表是亚当斯法。

数值方法的基本思想104104Euler公式1051Euler公式

欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点的一条称之为微分方程的积分曲线。积分曲线上每一点的切线的斜率等于函数在这点的值。105Euler法的求解过程是:从初始点P0(即点(x0,y0))出发,作积分曲线y=y(x)在P0点上切线(其斜率为),与x=x1直线相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线方程为当时,得

Euler公式106106同样,过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线的斜率=直线方程为Euler公式这样就获得了P1点的坐标。

107由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,…,Pn。对已求得点以=为斜率作直线

Euler公式当时,得108这样,从x0逐个算出对应的数值解Euler公式当时,得取109Euler公式从图形上看,就获得了一条近似于曲线y=y(x)的折线。通常取(常数),则Euler法的计算格式

i=0,1,…,n(2)还可用数值微分、数值积分法和泰勒展开法推Euler格式。以数值积分为例进行推导。将方程的两端在区间上积分得，110

选择不同的计算方法计算上式的积分项,就会得到不同的计算公式。(3)Euler公式用左矩形方法计算积分项

111

代入(3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式

由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。Euler公式1122梯形公式为了提高精度,对方程的两端在区间上积分得，改用梯形方法计算其积分项，即(4)

代入(4)式,并用近似代替式中即可得到梯形公式梯形公式113

(5)式的右端含有未知的yi+1,它是一个关于yi+1的函数方程,这类数值方法称为隐式方法。相反地,欧拉法是关于yi+1的一个直接的计算公式，这类数值方法称为显式方法。梯形公式

由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（5）比欧拉公式(2)的精度高一个数值方法。(5)

1143改进的欧拉公式显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。先用欧拉公式(2)求出一个初步的近似值,称为预测值,它的精度不高,再用梯形公式(5)对它改进的欧拉公式115改进的欧拉公式预测

校正（10）

可以证明,公式(10)的精度为二阶。这是一种一步显式格式,它可以表示为嵌套形式。校正一次,即迭代一次,求得yi+1,称为校正值,这种预测-校正方法称为改进的欧拉公式：116（11）或者表示成下列平均化形式（12）改进的欧拉公式1174

龙格-库塔方法1龙格-库塔(Runge-Kutta)法的基本思想Euler公式可改写成则yi+1的表达式y(xi+1)与的Taylor展开式的前两项完全相同,即局部截断误差为。改进的Euler公式又可改写成118

上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,而且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次f(x,y)的值,为一阶方法。改进欧拉公式需计算两次f(x,y)的值，它是二阶方法。它的局部截断误差为。龙格-库塔法的基本思想119

于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导,又提高了计算方法精度的阶数。或者说,在这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格—库塔（Runge-Kutta）法的基本思想。龙格-库塔法的基本思想1202二阶龙格—库塔法

在上取两点xi和,以该两点处的斜率值k1和k2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即式中:k1为xi点处的切线斜率值，k2为点处的切线斜率值,比照改进的欧拉 法,将视为，即可得二阶龙格-库塔法121二阶龙格-库塔法对常微分方程初值问题(1)式的解y=y(x),根据微分中值定理，存在点，使得式中

也即（13）K可看作是y=y(x)在区间上的平均斜率。所以122（14）将y(xi)在x=xi处进行二阶Taylor展开：

（15）二阶龙格-库塔法可得计算公式为：将在x=xi处进行一阶Taylor展开：123将以上结果代入（14）得：（16）

二阶龙格-库塔法对式(15)和(16)进行比较系数后可知,只要124式(17)中具有三个未知量,但只有两个方程,因而有无穷多解。若取,则p=1，这是无穷多解中的一个解，将以上所解的值代入式(14)并改写可得

二阶龙格-库塔法（17）

成立,格式(14)的