非线性振动漫谈

课程名称小组成员学院系别专业年级任课教师振动理论小组报告课程名称小组成员学院系别专业年级任课教师振动理论工程科学学院近代力学系理论与应用力学2012级陈海波2015年1月3日非线性振动漫谈1:概述牛顿方程建立以来，人们应用线性微分方程研究物体的运动，取得了许多重要的成果。实际上，不仅是力学，还包扌舌其他物理学领域，人们长期研究的主要是线性理论。自20世纪60年代以来，人们在实验研究中发现，在自然科学的许多领域都存在一些线性理论无法解释的现象，开始是非线性振动，后来又在流体力学和声学领域，接着是伴随激光而诞生的非线性光学...非线性向人们展示了一副令人惊奇、甚至是不可思议的图景。本文首先宏观得介绍非线性现象及其与线性系统的不同点，然后初步介绍非线性振动的两类求解方法，一类是定性的方法或称几何法，用以判定一个系统的发展趋势，主要是其稳定性问题；另一类是定量的方法，主要有平均法，迭代法与摄动法等。由于定量法设计的数学知识比较复杂，且数学推导很麻烦，因此本文将把重点放在定性求解上。2：非线性系统与线性系统的对比从数学的角度看，线性系统有两个显著的特点：一个特点是，因为自变量与函数之间的关系是线性的，因而自变量的变化率与函数的变化率之间成确定的比例。例如，函数夕=ax+b,它对自变量x是线性的，则由可见，当时，也有另TO。这意味着函数值对自变量的取值精确度不敏感，亦即相应于自变量的微小变化，函数值也只会产生微小的变化。而一般的非线性系统，对初始值具有高度的敏感性，著名的蝴蝶效应就是典型的例子,意思是说，今天在巴西的热带雨林里有一只蝴蝶偶尔拍打了一下翅膀，可能在两个星期后,就在美国的德克萨斯州引起一场龙卷风。天气的演化对初始值如此敏感，因此长期的天气预报是不可能的。线性系统的另一个特点是，系统的整体性质可以由组成它的各个子系统的代数叠加得出，这就是所谓的线性叠加原理。如杲一个系统可以分解为若干相对独立的子系，则线性关系表明的是，只要各个子系的行为都已知，则系统的所有行为就是所有子系的简单叠加。从这里可以看出，线性系统是由若干互不相干的独立子系组成的，线性关系就是来源于各个独立子系的独立贡献。于此对照，非线性系统的各个子系之间有着不可忽略的相互作用，因而在非线性的情况下，爭情就变得复杂多了。下面以洛伦兹方程为例子，简单得说明非线性现彖对初值的高度敏感。1963年美国麻省理工学院的气彖学专家洛伦兹在研究天气预报问题时，将异常复杂的人气对流偏微分方程化简后，得到一组常微分方程X=—10x4-1Oy•(1)<y=28x—y—xz•Sz=-_Z4-xy这一方程已成为混沌理论的经典方程，即著名的洛伦兹方程。洛伦兹当时使用的是真空管电子计算机，通过数值计算求解这一方程，因为计算费时太多，为避免每次从头算起，他把每次计算的中间结果打印下来。使他惊异的是，在重复上次计算结果的过程中，知识开头一小段与上次结果一致，但很快就越来越偏离上次的结果。他意识到，他的这一组方程并不是线性的故不遵从传统线性理论下的规律，而是对初值具有高度的敏感性。除了对初值极为敏感这一点以外，洛伦兹方程的解(图2)也显出很奇异的特性。从图2可以看到，方程的解总体由两个坏套组成，像三维空间中的某种双螺旋，又像蝴蝶的两个翅膀，这种结构称为奇怪吸引子，因为它们吸引相平面的相点，但所有的相点又都永远到不了坏套的中心。从图中看，每一个环套上都分布着细密的轨线，轨线一层层缠绕，在一个坏套上转几圈，又跑到另一个换套上，完全无法预料什么时候从一个坏套跑到另一个环套上。而且,这些轨线尽管紧密缠绕，但从不相交。从轨线的这些特征可以明显看到，洛伦兹方程的长期行为是无法预测的。图2：洛伦兹吸引子3：非线性振动的求解方法3.1非线性振动的定性分析方法3.1.1状态空间及状态方程将广义坐标q（t）与J'•义速度q（t）组成一个向量{q1,q2,...,qix; »称为状态向量,其各个分量称为状态变量，状态向量所存在的空间，称为相空间或状态空间。状态向量的端点，称为“状态点”，其运动轨迹称为相轨迹或轨线，轨线的总体称为相图。通过相空间中一点，一般只会有一条轨线，因此其各条轨线不会相交（个边点除外），而整个相图纹理井然，便于分析，因此非线性系统的几何理论经常在相空间展开。以状态变量作为基本变量可得到状态方程，多自由度的运动方程一般可以写为q(t)二^(q1q2,q3,...,qn；qpq2,q3,-,qn；t),i=l,2,3,...,n为此令xi=Qi>耳+乂=Qi,i=1,2,…,n（3）而状态向量成为(4)因此多自由度的运动方程一般可以写为Xx（t）=普（耳卷,…,Xk；t},i=l,2，…⑴（5）方程（5）称为系统的状态方程，这是2n个变量的一阶微分方程组。对于一组确定的初始条件

x1(0)=«1,i=l,2,...,2n(6)可以证明方程(5)有惟一的解。这表明通过状态空间中一个确定的点，一般只有一条确定的轨线。3.1.2平衡点及平衡点的稳定性状态空间中{x}h{°}的点称为普通点或正则点，而••••<x}={Xl,X2,...,X2n}T={Xi(^,X2>...,X2n)}={0}⑺的点称为奇点或平衡点。在平衡点上，所有的状态变量的变化率対(1=1,2…，2门)均为0,因而状态变屋不会改变。其结果是系统只能静止在原來的位置上，不可能运动。如果平衡点不在原点，那么按下式进行坐标平移：(8)就可以将原点平移到平衡点上。平衡点可以分两类,即稳定的平衡点和不稳定的平衡点，其差别并不在于平衡点本身的状态,而在于系统在略为偏离平衡点时的运动趋势是趋向于回到平衡点，保持在平衡点附近运动还是趋向于偏离该平衡点越来越远。相应的，该平衡点分为渐进稳定，仅稳定的或不稳定的,其中前两种平衡点又称为稳定的平衡点。3.13单自由度自治系统状态方程在平衡点领域的线性化。如果状态方程(5)的右边不显含时间t,即Xi(t)= = (9)则状态空间的流场是定常的，这种系统叫做自治系统。单自由度系统的状态方程式如式(10)所示(10)X1=Xiaq,X2),X2=乂2(络，冯)(10)设平衡点的坐标为X1=apXj=a2代入式(10)得(11)Xi= a?)=0,X2= (a】，a?)=0(11)由于XI,X2是非线性函数，一般会得到关于31/32的多组解，即可能有多个平衡点，我们总可以采用(8)所示的坐标变换，将坐标原点移到任何一个平衡点，因而不失一般性我们可以认为平衡点即为原点，即al=a2=0,而在原点附近将状态方程式(10)展开为式中(12)ay=式中(12)ay=(13)勺是二阶以上的微量。如略去这些高阶微量，并采用矩阵记法有

此式在原点(平衡点)附近近似地成立。一下分析此式所代表的近似线性系统在原点附近的相图的几何特性，并按之对平衡点进行分类。3.1.4平衡点领域中的相图及平衡点的类型式(14)表明系统在平衡点附近的动态特性由矩阵⑻确定，这里假定[a]是非奇异矩阵。为了研究次矩阵对平衡点附近的相图性质的影响，对相平面进行线性变换，其目的是将⑻变为尽可能简单的形式。若记{x}二[b]{u}(15)其中[b]是一非奇异的变换矩阵：｛u｝是新的状态变量。对于新的状态变量｛u｝的状态方程[a][b]{u}=[c]{u}(16)[a][b]{u}=[c]{u}(16)其中[c卜[b「】[a][b][c],[a]两矩阵称为相似矩阵。相似矩阵具有相同的特征值。矩阵[a]的特征值满足以卜方程:a】i_兄an二。a21玄三一兄展开上式，得到关于入的二次代数方程组：尤-(尙+a22)；l+(aiia22_ai2a2i)=0(17)求解可得出两个特征值入1，入2.这两个特征值有以下三种情形，对于每一种情形，都可以选择适当的变换矩阵[b],使得[c]矩阵化为最简单的约当(Jordan)形入1,入2是相异的实数，此时有[c]=[b「[a][b卜卜.(18)式(16)成为

其解为1片=1加工=1屛（20）如果入2〈入1〈0则在原点附近的相图如图2所示。这种情况下，相平面（ul,u2）的原点称为节点。从相图上可见，从原点附近的所有点出发的轨线都单调地趋向原点。因此，此节点是渐进稳定的。如果入2>入1>0,轨线的形状任与图2相同,但所有的箭头需反向，我们得到不稳定的结点。如果入2<0<X1则相图如图3所示，这时的原点称为鞍点，鞍点总是不稳定的。图3鞍点（2）如果入1,入2是相等的实数。则c的约当标准型有两种:[c]=AA[c]=AA对于第一种情况，解为4=4001112=1/6创(22)这时原点仍为结点，但相轨迹成为过原点的直线，这种结点称为边界结点。且当入1<0时是稳定的，而当入1>0是不稳定的。至于式(21)所示的第二种情况，将得到一种退化的结点，轨线为曲线，仍然是入1<0时是稳定的，而入1>0是不稳定的如果入1,入2是共轨复数，令入1二(a+ip)，入2二(a-ip),a,0为实数，其解为做线性变换做线性变换以(23)式中的ul,u2代入式(24),并记ulOu2O=uO,得^(iioeTsinQt'V^QiocTcosQt⑵)在(vl,V2)平面上这时一对数螺旋线，如图4所示，0的符号确定螺旋线的旋向：0>0,为逆时针方向；P<0,为顺时针方向。a的符号则决定是向内旋，还是向外旋，即决定平衡点的稳定性：a<0,向内旋，渐进稳定；a>0向外旋，不稳定。图3是a<0,0<0的情况，这种类型的平衡点称为螺旋极点或焦点。在平衡点为焦点的情况卞，其附近的轨线以衰减振荡的方式趋于平衡，或以增幅振荡的方式，偏离平衡点。当a二0时，轨线称为图5所示的同心圆，这种平衡点称为中心，属于仅稳定。图5中心图图5中心基于以上的讨论与分析，可以不必求出［a］的特征值，而直接按［a］的元素来判断平衡点的类型。若记坷|+%2=p,ana22-a12a21=q,A=p2-4q奇点的不同类型由参数P和D完全确定，只要这两个参数确定了，则系统奇点的类型就确定奇点类型和这两个参数的关系可以归纳如卜:A>0q>0结点<A<Q<q<0A>0q>0结点<A<Q<q<0鞍点p=0中心pHO焦点“p<0p>0p<0稳定结点不稳定结点稳定焦点不稳定焦点3.2非线性振动的定量分析方法3.2.1解析求解法非线性振动问题很难严格求解，但有限摆幅的单摆问题却是一个成功的例子。设单摆的质屋m,摆长为I。选取摆线偏离铅垂线的角度e为广义坐标，以8=0点为势能零点，体系的拉格朗口函数为L=imZ262-mgl(l-cos0)=imZ2e2-2mglsm2|因为器=0,所以广义能量枳分（这里就是机械能）守恒，即otimZ262+2mglsm2-=E022设初始条件为e=e°且e=o,代入上式得Eo=2mglsin2—2因而di为了求得系统的周期T,我们将上述方程在te[o£]和66[00/0]为了求得系统的周期T,取负号，可变形为•ode0o•ode0odugJodugJo心一以sin%全亍K(k)(2)d9=(1)e2做变换sin|=Ksinu,06[0,0O],uG[o冷]其中K=sil占则式（1）变为其中z、f2duK(k)= 「JoVl—K2sin2u是第一类椭圆积分。式(2)表明，人振幅单摆的周期与振幅有关。为具体考察T与6。的关系，将式(2)展开成以的幕级数,得KX)=畀+(霁2+(拎2/+£.持+…]所以T=2ir匸(1+-sin2—+—sin4—H—)y]g\4 2 64 2 )可见，当振幅60«1时，上式中的S曲2単以后的高次项均可略去，得到与振幅无关的小振动的周期T=2ir然而随着振幅的增人，sin2^以后的高次项不能忽略，T随％的増人而增大。下表是具体的数值计算结果。表1：不同摆幅下严格解与零阶近似的周期比较0o/°1257.22422.81304590 180T/1.00001.00001.00041.00101.01001.01741.03991.1803 8T。28800174由上表中数据可以看出，在精度要求一般时，比如1%精度，要求初始摆动角不超过22.81°即可，这是一个很容易满足的要求。但如果要求精度再增加一个数量及，如在用单摆测量重力加速度实验中，则％必须小于7°,实际操作时务必注意这一限制。此外，如呆取％取180°时，表中的周期为无穷大，但这仅有数学上的意义。因为它对应一个不稳定平衡的初始位置，而实际情况下总会有各种因素打破这一不稳定平衡。3.2.2.微扰法微扰法又称逐级近似法摄动法，它适于求解形如X+co02x=f(x,x)的振动方程，其中非线性项|/(%,%)|«td02%,可将其视为对线性方程的扰动或摄动，为表示微扰的数量级，我们引入参量£,将上式写成x+cu02x=ef(x,x)(3)微扰法的基本做法是：令%和s的n阶近似解为(%=Xo+ +E2X2+■■•+EnXn(co=O>0+£(jOl+£2(j02H F£nd>n其中£%和0©(i=l,2,…,n)为对零阶近似解勿和如的i阶微扰。将上式代入式(3)中,使等式两侧£同级幕的系数相等，就得到关于各阶％和3的一组方程，从零阶开始逐级求解这些方程(求解第i阶方程时，把零到i-1阶的解作为己知条件)，便可以得到原方程的n阶近似解。下面我们仍以人摆幅摆动的单摆为例，说明微扰法的应用。从单摆的拉格朗口方程易知其运动微分方程为d+ysill0=0 (4)将sine做幕级数展开sin&=e-抨+盒沪一…取到e的三次项，式(4)近似为d+a>o26=a03其中如=占,a=ajQ2/6.上式可以看成是线性方程0+too20=O的解再加上小量a/的微小修正后得到的。为了明确表示数量级，引入£，将上式变成O+a)Q2O=Ea03(5)假定只求二阶近似解，则令0=Oq+£&]+8^02co=a)0+ecoi+e2co2代入式(5),得@()+801 — —£2£()2)2(&0+£0^+异&+ &2)'略去£的三次以上项，比较上式等号两侧同为日，』和E2的项，分别得到零阶，一阶和二阶(歹0+CO20q=0近似方程分别为｛ +co20l=2COCOX0Q+aOQ3 (6)VO2+ =-CDi^Oq+2C0C01&1+26O6O2&0+3cl0q^0零阶近似即常规的小振动单摆的结果为00—力0cos(a)t+卩)其中4。和。为待定常量。将上式代入式(6)中第二式，则01+a)20i=2a)a)iAQcos(a)t+卩)+aA03cos3^a)t+<p)利用三倍角公式化01+0)28]=3 3(2coco1j4o+乂:J)cos(o)t+<p)+弓-cos3(et+<p)(7)如果我们将该方程的右侧部分视为策动力项，则其中的第一项的频率与左侧的固有频率都是3,由于这里不存在任何耗散，此策动力将产生振幅无穷人的共振，除非该项的系数为零,即3aA032&)0)^0H =0所以3aAQ23=_ 丄 8a)于是式(7)简化为.. -ciAq^0i+a)2址=—-—cos3(cot+(p)4易解其特解为(通解己经包含在零阶近似解中，不必再考虑)3%=AiCOs3(3t+0),Al=-^采用类似但繁琐的步骤，可以得到二阶近似卞的求解结果。首先，为去掉不符合实际情况的发散型共振，要求—3ciAqAy2Aq^8a)8a)因而二阶微分方程简化为。2+ =“331+严0Ticos3(3t+y)+3fl^ cos5(o)t+(p)解得02=i42cos3(cot+<p)+F2cos5(cot+<p)A_(4cucui+3a4o2)4iD_aAo2Ai&2=02=_-系数可进一步简化为=3o%ip二a%52 — 1024- 1024』于是我们得到二阶近似下的总解为8=&0+&1+&2=Aqcos(cot+cp)+(4丄+i42)cos3(cot+<p)+F2cos5(cot+<p)=久cos^t+0+(-；；：；l：；/：：)cos33+初+需咅cos5(加+0)一般情况下没有必要求解更高阶近似，因为作为问题的出发点的运动微分方程中仅精确到a的一阶项。最后，有初始摆角为弘和角速度为零的条件，即卜谿cos3叶議討S50"。Aqsin</>+(-sin5°=03ai4028co3a4o3,9a24o\・o,5a240532w2+1024cu4^Aqsin</>+(-sin5°=03ai4028co以及频率关系21a%4-4-256to3联合求解各阶振幅，初位相和频率（实际计算时町用迭代法求快捷解）。卜表是不同初始摆角下二阶近似解与解析解周期的比较。由表可知，直到初始摆角达到人约15°,二阶近似解与严格解的误差才有10-5,即使摆角增加到90°,误差也仅为1・6%°表2：不同摆幅下二阶近似解与严格解的周期比较00/°14.2520304590T1.000011.000041.000201.001031.01576/%值得指出的是，微扰法中引入£来表达阶数，只是为了方便区分各阶小量，所以一旦各阶量求得后，在表达总量时就可以舍弃它了。3.2.3求解非线性方程的数值分析方法数值分析方法就是对非线性系统进行数值积分，在时域内把响应的时间历程离散化，对每一时间步长内可按线性系统来进行计算，并对每一步的结果进行修正。这种方法又称为逐步积分法或直接积分法。数值分析方法得到广泛应用一个原因是因为非线性分析理论发展的不完善性，对很多问题无法进行理论上的分析；另一个原因是数值分析理论的发展和计算工具性能的提高使得数值分析成为可能。常用的数值分析方法有纽马克（Newmark）法、威尔逊（Wilson）法、Runge—Kutta法等。由于篇幅有限这里不再详细介绍。3