动力学问题的有限元法.pdf

1、第七章 动力学问题的有限元法 第七章 动力学问题的有限元法 结构动力学是研究动载荷作用下结构动力反应规律的学科，讨论结构在动力荷载作用下反应的分析方法，寻找结构固有动力特性、动力荷载和结构反应三者间的相互关系。研究结构在动力荷载作用下的反应规律，能够为结构的动力可靠性（安全、舒适）设计提供依据。前面介绍的静力学问题的研究对象是受不随时间变化的载荷作用。而动力学问题的对象受随时间而变的载荷的作用，从而使在结构中产生的位移、速度、应力和应变都随时间而变。当结构受随时间变化的载荷作用，且这种载荷的作用对结构的变形和应力的产生起主要作用，以致影响设备的安全性，或舒适性。这时就要进行动力学分析，充分认识

2、其规律性，从设计阶段就抑制这种不利状况的发生。例如，有时虽然动载荷不大，但结构在交变力的作用下，其某些固有频率与激励力的作用频率相接近时，就会引起很大的振动、变形或应力，这时，就必须对结构作动力学分析。又如，要利用结构在周期性作用力驱动下的定向振动，例如利用这种运动输送产品，这时，就必须巧妙地设计结构，使其具有某些与激励频率一致的固有频率，并且使结构对激励具有适当的响应能力。总之，不管是利用振动，还是抑制振动，都需要进行结构动力学分析。当前结构动力学的研究内容有三类。第一类问题：反应分析（结构动力计算），第二类问题：参数（或称系统）识别，第三类问题：荷载识别。第一类问题是已知系统动态特性和动载

3、荷作用部位及大小，求出系统的响应随时间变化的位移，速度，加速度和应力等。第二类问题是已知系统的输入输出特性，分析系统固有的动态特性，结构模态分析就属于这一类问题。第三类问题是在已知系统动态特性的条件下,通过测量系统的响应,或由响应准则预先给出响应要求,以此识别对响应的外载荷。三类结构动力学研究内容的载荷、结构和响应之间的关系如图 7-1 所示。动载荷种类大致分类如图 7-2 所示。图7-1 结构动力学研究的内容 图 7-2 动载荷种类 1本章主要介绍结构动力学分析的基础知识，并主要介绍系统固有特性的有限元分析方法有限元模态分析。主要知识点：1.导出有限元的结构动力学方程和模态分析有限元方程。2

4、.介绍模态分析方程中质量矩阵与阻尼矩阵的计算。3.把计算结构的固有频率与振型的问题归纳为一个求特征值的问题。7.1结构动力学方程及有限元方程 7.1结构动力学方程及有限元方程 系统的受迫振动的微分方程：m()xcxkxR t+=?其中：质量，阻尼系数，弹性系数。mck,x x x?分别为加速度，速度和位移。1.单元的动力方程：（用“虚功原理”求解）结构动力学有限元分析与静力学的有限元分析一样，首先要建立单元的有限元方程。为求解起见，要引入达朗贝尔原理。质点的达朗贝尔原理可表述为：当非自由质点运动时，主动力，约束反力和惯性力构成一个动平衡力系。惯性力是一种虚构的力，由质点的加速度引起。mFNS

5、即：0fns+=现考察单元的动力学方程的建立。1).准备工作：将位移函数 表达成近似函数。单元内任一点位移 (,)()eNN x y z=et，是单元内位移的近似函数,这儿 N只是位置的函数，与时间无关。而 e与时间相关。(,)()eeNN x y zt=eeNN=?2).分析单元上的作用力。图 7-3 单元分析 主动力：集中力，面力，体积力：psRRRRV=+阻尼力：与速度相关，并与速度方向相反，两者的比例关系为阻尼系数：ceRCC N=?2惯性力：与加速度，质量相关：单位体积惯性力 meRN=?,阻尼力、惯性力与主动力方向相反。密度，单位体积质量。3).建立动力学平衡方程：利用虚功原理和达

6、朗贝尔原理，系统处于动平衡，首先将载荷移置到节点上去。利用虚功原理。外力做的虚功：11*11*TeTTppTeTTSsSsTeTTVVVVTeTTVCVCTeTTVmVmRNRRdSNRdSWRdVNRdVNRdVRdVNRdV=集中力：面力：体积力：阻尼力：惯性力：RdV 内力做的虚功：*TeTTeVVUdVBDB=dV0 虚功原理：()WUWU=，最小势能原理 根据虚功原理，应有 WU=1*1*eTeTeTTTpSsVVeTeTeTTTTeVCVmVNRNRdSNRdVNRdVNRdVBDBdV+=由*eT的任意性，可约去，且由阻尼力，惯性力的定义可知：cmeeRCC NRN=?，代入上式

7、，有 11TTTpSsVVeeTTTeVVVNRNRdSNRdVBDBdVC NNdVNNdV+=+?=eRi令：11TTTepSsVVTeVTeVTeVNRNRdSNRdVRKBDB dVcC NN dVmNN dV+=上式可表达为：，这就是单元的动力学有限元方程。eeemcK?+=+=2.结构整体动力学有限元方程 与静力学的集成原理相同：根据 局部单元节点编号与整体的节点编号的对应关系，利用，即一个节点上的节点力是该节点所在的所有单元的相应节点（同这个节点）上的节点力的总和来集成。eeiiRR=3位移唯一性。得 MCK+=?R eMm=，eCc=eKk=其中，当 0C=，称为无阻尼受迫振动

8、，此时 MKR+=?当 0,0CR=，称为无阻尼自由振动，此时 0MK+=?7.2 单元质量矩阵和阻尼矩阵的表达式7.2 单元质量矩阵和阻尼矩阵的表达式 在动力学有限元方程中除刚度矩阵之外，还有两个重要的系数矩阵：即质量矩阵和阻尼矩阵。这一节中将介绍单元质量矩阵和单元阻尼矩阵的表达式。一.单元质量矩阵 单元质量矩阵有二种表达形式：一致质量矩阵和集中（堆积）质量矩阵。1.一致质量矩阵：用公式 计算的单元质量矩阵，称为一致质量矩阵，是因其中的 TeVmNN=dVN就是位移模式形函数 N而得名。或者说，推导质量矩阵采用的形函数与采用推导单元刚度矩阵时采用的形函数是一致的。以平面应力问题为例：三角单元

9、，ijmNNININI=m（位移模式中的 0 0 00 0 0 ijmijNNNNNNN=N）图 7.3 三节点三角形单元 则三角形单元的一致质量矩阵 em：（以单位厚度计）ii22jj2j 0 N 0 N 0 N 0 0 N 0 0 N 0 0 jmijmijmeTijNNNNNNNNmN NdxdyNN=?j2mm2mm NN 0 N 0 0 0 N 0 N 0 mijmijmdxdyNNNNNNN?由于m，利用数值积分：2 iNii 0 0 N 0 N,iijjmNL NLNL=4!2(2)abcijma b cL L L dxdyabc=!+?，有22!24!62(1 12)!12ii

10、jN dxdyN N dxdy=+从而得到 20 2 1 0 212 0 1 0 21 0 1 0 20 1 0 1 0 2em=对称 2.集中质量矩阵：将单元质量假想地集中到节点上去，也就是说只在节点上有质量。这样，一个节点的加速度就不会影响其它节点的惯性力（相互之间不受影响）具体处理方法是将质量平均分到所有节点上去。例如三角形三个节点单元，每个节点分到 1/3的质量，即1m3it=,t为厚度，为三角形面积。11 1 1 1 11 1 133 1 1 1 1ewmt=资料表明，集中质量矩阵计算出来的频率与一致质量矩阵的计算结果相差无几。在单元相同的情况下，集中质量矩阵计算的结果低于一致质量矩

11、阵的计算频率，而振型则是一致质量矩阵的计算结果精确。由于集中质量矩阵式对角线矩阵，所以计算中常采用。3.附加质量 当有附加的质量在节点上时，要在计算整体质量矩阵时，加上这些附加质量cm 即 ecMmm=+二.单元阻尼矩阵 根据阻尼的成因可以分为粘滞阻尼和结构阻尼。1 粘滞阻尼（比例阻尼）前面得到的阻尼矩阵 是将阻尼看作是正比于节点的运动速度，称为粘滞阻尼。它消耗振动的动能（能量）。从 可知，这个阻尼矩阵与质量矩阵成正比。TeCC NN d=VdVeTmN N=eeCCm=，/Cc=如果,c均为常数的话。所以称为比例阻尼。（相对于质量矩阵而言）NTeCC Nd=V2结构阻尼 比例于结构的变形速度

12、?，由于材料内部的摩擦引起的，称为结构阻尼。这时的阻尼力可以简 5化为 cRD=?，这样可得到单元的结构阻尼矩阵使其具有 的形式，所以结构阻尼正比于刚度矩阵：BeTkCCBDd=V eekCk=3Rayleigh 阻尼 在实际分析中，精确决定阻尼矩阵是相当困难的，通常允许将整体阻尼矩阵简化为质量矩阵M和刚度矩阵 K的线性组合，这样可用来解耦有阻尼的振动方程。M+KC=，称为 Rayleigh 阻尼。其中j22i222(-)2(-)ijiijjijjiji =其中22iii+=称为结构的第i阶振型的振型阻尼比振型阻尼比。2i+称为结构的第i个校准振型的振型阻尼常数。工程中一般取i=0.020.2

13、4，它与结构的类型，材料的能量消耗特性及振型有关。可通过实测，从系统的固有频率和辐频特性曲线获得。质疑质疑：相关内容的矛盾描述：在王勖成的书 P476 有如下字样：“,是不依赖于频率的常数”。但是，在同一本书的 P487，又出现下列字样：“如果根据实验或相近似结构的资料，已知二个振型的阻尼比i，j，可以得到 j22i222(-)2(-)iiiijjijjiji =”7.3结构动力学特性的计算7.3结构动力学特性的计算 求解结构的固有频率和对应于每一阶固有频率的振型，这是动力分析的基本内容，称为模态分析。振型是指在某一个自振（固有）频率下，结构各节点振动值的集合，即节点振幅的相对大小。由于阻尼对

14、结构的自振频率和振型的影响不大，所以，一般只考虑无阻尼的自由振动。一.模态分析的有限元方程 在 MPCK+=?中，令 P0,0C=，得模态分析的有限元方程：6 M0K+=?（注：注：联系到温度场的瞬态热传导方程：3PKtKt+=?，P0=无内热源。有人作比对分析，对温度场的分析采用所谓的“热模态分析”，但要注意，两个方程只是“形似”，它们的解有区别，热是一阶方程，振动是二阶方程（相对于时间）热是指数型解，振动是正弦函数分布。）二。数学预备知识 二。数学预备知识 求解固有频率和振型，数学上涉及到矩阵的特征值问题。这儿做一个有关矩阵的特征值，特征向量和特征方程方面的知识的扼要回顾，以便将我们的主要

15、精力放到模态分析的物理意义上去。1 特征值和特征向量 对 n 阶方程：和 n 维非零列向量，如果有一个数1111nnnnn naaAaa=?12(,)Tn=?使得矩阵n nA和维非零列向量n有A=，则称为矩阵A的特征值（特征根），向量为矩阵A的特征值所对应的特征向量。()0AAI=，为非零列向量，若等式成立，应有0AI=2.特征矩阵 11121212211 nnnnaaaaaaAIaan=?称为特征矩阵，其中I为 n 阶单位矩阵 3.特征多项式：行列式AI称为矩阵A的特征多项式。记为()AI=4.特征方程式：()0=，称为A的特征方程。特征方程()0=的个根n12,n?就是矩阵A的个特征值。n

16、集合12,n?称为矩阵A的谱，记作 chA线性齐次方程组()iAI0=的非零解便是矩阵A的特征值i所对应的特征矢量。在线性齐次方程组()AI0=中：AI特征矩阵 i特征值 特征向量 0AI=特征方程 三.自由振动的特征值和特征向量 自由振动方程 MK+?0=（1）1.求特征值和特征向量：7自由振动 各点位移为平衡位置附近的简谐振动：0cos()wt=+这样，20cos()t2=+=?，其中 0为振幅，为频率（谐振频率），为相位移。将它们代入方程（1）得：20(M)K+=0 （2）由于，即振幅不可能全为 0，为使式（2）成立，必须有 00 2M0K+=即，要使 0得到非零解，必须要有 2M0K+

17、=，从中可以求出i（i=1，n）,n 为自由度数。然后，将i代回到式（2）求出每个i对应的 0。(一般工程上只需少数个最低的谐振频率即可)，对应于每个ni，都有一组 0i 记 2ii=广义特征值工程上称为谐振频率，自振频率。0i特征向量工程上称为结构的振型，或振动模态。2 特征向量的性质 1）0解不是唯一的，当 0为解时，0 也是解。将 0 代入（2）式时，有 20(M)K+=0 0 也是结构的一组特征向量，也就是说，在某个频率下，物体的各点的振动相互之间保持一个相对位移的比例关系，但它们的绝对值是可变的。2).0的归一化处理 通常将 0归一化处理，取 0使得 001TiiM=（3-1）这样归

18、一化的 0具有以下性质 000010 TiiTjiMM=（3-2）0000 0 TiiiTjiKK=（3-3）证明：由 001TiiM=，有其余的表达式成立。8 证明有：00TiiiK=证明过程:由有限元方程 0()iiKM=0 有 00iiiKM=左乘得 0Ti 00000()TTTiiiiiiiKMM=0i 001TiiM=001=TiiiiK=证明有：及 000 TjiM=当 000TjiK=ji 证明过程：1）.对 阶i 0iiiKM0=a)对j阶 00jjjKM=b)对 a)式左乘 0Tj 有 000TT0jiijKiM=（1）对 b)式左乘 有 0Ti 000TTijjiK0jM=

19、（2）2）.将（1）式转置：得 其左边：00000()TTTTT0jiijiKKjK=（3）其右边：000()TTTijiiiM0jM=（4）这儿利用了 K和M的对称性：，TKK=TMM=（3）=（4）有 000TTijiiK0jM=（5）3）.将（2）-（5）有：000()TjiijM=当ij时，有 000TijM=代入（5）式 有，得证。得证。000TijK=3 振型矩阵（模态矩阵）和频谱矩阵（特征矩阵）将 001,n?集成X方阵，n nX称为振型矩阵或模态矩阵。9这时 1001Tn nXMXI=?1200TXKX=?，特征矩阵，频谱矩阵 X是关于 ,MK的正交矩阵。例一：如图示悬臂梁，不

20、计轴向位移，求其固有频率和振型。图 7-4 悬臂梁 1）位移模式 无阻尼自由振动方程：0MK+?=,模态有限元方程：()200KM=特征方程为：20KM=；不计纵向（轴向）位移，有：ieijjjviv=，其中挠度,转角21234vaa xa xa x=+3dvdx=边界条件：110v=，所以 223300evv=2）单元质量矩阵 取两个单元。单元质量矩阵 eTVMNN dV=232323232322321 32232xxxxxxxNxllllllll=+x 10 23232323232323223223223232321 3221 3223232elxxllxxxxxxxxxxxllMxAll

21、llllllxxllxxll+=+dl+考虑到一维单元的自然坐标 1ijxxLLLLL=，将eM转换成 L表达的形式，可以使积分变得方便。2322322322321 32(12)1 32(12)32()32()jjejjjjjjjjjjjjljjLLxLLMLLxLLLLxLLAdLLxLL+=+l+11112121111121211212222212122222xxyxxxyexxyyyxyyyxxyxxxyxxyyyxyyMMMMMMMMMMMMMMMMM=y 23211(1 32)xxjjlMALLdl=+=1335lA()23221111(1 32)(12)210yxjjjjlMALL

22、 xLLdllA=+=()232321(1 32)(32)xxjjjjlMALLLLdl=+=970lA()23221(1 32)()yxjjjjlMALLxLLdl=+=213420lA()2211(12)yyjjlMAxLLdl=+=21105lA 23222(32)xxjjlMALLdl=1335lA()2222()yyjjlMAxLLdl=+=3105lA()()232121332(12)420yxjjjjl2MALLxLLdlAl=+=()32212()(12)140yyjjjjllMAxLLxLLdlA=+=()()2232211()32210 xyjjjjl2MAxLLLLdlA

23、l=+=1122311112121222111121211212222222121222223322131191335210704201111321010542014091313117042035210131142014021010 xxyxxxyxexyyyxyyyxxyxxxyxxyyyxyyyllllMMMMllllMMMMMAMMMMllllMMMMllll=5 2222156225413224133541315622420133224elllllllAMllllll=3）单元刚度矩阵 梁的刚度矩阵可以用两种方法求得。利用单元梁的平衡，利用能量法。利用单元梁的平衡。这儿要重温材料力学中

24、关于梁的理论：纯弯梁微分单元受弯矩M和剪力V的作用，梁的微分单元平衡时，有 2323,=d vd vMEIV EIdxdx=。图 7-5 梁的微分单元 图 7-6 梁的有限元单元 从图 7-5 和图 7-6 的比较中可以看出，,ijMMMM=；,ijVVVV=单元位移模式：231234evaa xa xa xN=+=232323232322321 32232xxxxxxxNxllllllll=+x ieeeijjMVkRMV=，利用公式232,=d vd vMEIV EIdxdx=3，和位移模式 evN=12可得 33333232126126eeidNd vVEIEIEIdxdxllll=22

25、222322321261266462eeiiiiidNxxxxd vMEIEIEIlldxdxllllll=+33333232126126eeidNd vVEIEIEIdxdxllll=22222322321261266462eejjjjjxxxxdNd vMEIEIEIlldxdxllllll=+现令 0,ijxxl=，有 323222223323222221261261261266462646212612612612662646264iiieeeijjjjjjllllvvVlliiMllllEIllllkREIvvVlllMllllllllllll=22322126126646212612

26、66264ellllllEIKlllllll=，式中，l为单元长度，为单位长度质量。利用能量法求单元刚度矩阵 TeVKBD B dV=，几何矩阵 B由几何方程求得。,eduBdx=022bdudxd vydx=是拉伸引起的正应变0dudx=和由弯曲引起的正应变22bydydx=v的总和。本题不计轴向位移。00dudx=2222ebdNd vyydxdx=22dNBydx=由 evN=和 232323232322321 32232xxxxxxxNxllllllll=+x 有 22232236124661226dN2xxxxxByylldxllllll=+13 232223223223261246

27、612466122661226TeVVxllxxxxxllKBDB dVyEdVxllllllllllxll+=+22323223232322223222322612612466126126122646612464661246eVxxxxxxxllllllllllllllxxxxxxllllllllllllKEy +=2223232322323222322222661261261246612612262661226462661xllxxxxxxxllllllllllllllxxxxxlllllllllll+232226dVxxlll+223232232323222232223226126124

28、66126126122646612464661246exxxxxxxllllllllllllllxxxxxxllllllllllllKEy dA +=2223232322323222322222661261261246612612262661226462661xllxxxxxxxllllllllllllllxxxxxlllllllllll+232226Dldlxxlll+其中 抗弯刚度。2DEy dAEI=令 2232322323232222322232612612466126126122646612464661246xxxxxxxllllllllllllllxxxxxxlllllllllll

29、lF +=222323232232322232222326612612612466126122626612264626612xllxxxxxxxllllllllllllllxxxxxxllllllllllll+2226ldlxll+对其中的每一个元素分别积分，得到：2112361212lxFd3lll=+=l，122323612466lxxlFdlllll=+=l，132323361261212lxxFdllllll=+=14232612266lxxFd2lllll=+=l，2222464lxFdllll=+=，232233466126lxxlFdllllll=+=14142246262lxx

30、Fdllllll=+=，23323361212lxFdlll=l，342322612266lxxFdllllll=+=2442264lxFlll=+=dl 其他各元素可利用对称性求得。2232212612664621261266264ellllllEIKlllllll=4）质量矩阵和刚度矩阵的集成 Mmm=+，Kkk=+222223221261260064620012624012662086200126126006264llllllllEIKllllllllllll=，222222215622541300224133005413312054134201330813300541315622001

31、33224lllllllllAMlllllllllll=（注：注：当有较多位移边界条件已知时，可以在集成计算各元素之前，先作降阶处理，划行划列，划去相应的行与列，只集成剩余的各元素。这样可节约大量时间。）5）求固有频率 代入边界条件：110v=，所以，22223120541308133541315622420133224lllllAMllllll=，2232224012608621261266264llllEIKlllllll=特征方程 20KM=22222322222401263120541308620813301261265413156224206264133224llllllllEIlA

32、lllllllllllll=令：42420lAEI=则方程化为222224312 0 (1254)(613)0 (88)(6 13)(23)0(1254)(6 13)12 156(622)(6 13)(23)(622)(44)llllllllll+=+得：()443217880114298481170936803521440l+=1512340.001841,0.073481,0.840564,7.080982=求得：143.52(2)EIAl=，2422.22(2)EIAl=3475.16(2)EIAl=，44218.14(2)EIAl=1223.523.52(2)EIEIAAlL=,222

33、2.22EIAL=,3275.16EIAL=,42218.14EIAL=这儿。2lL=而精确计算值为 143.515(2)EIAl=，2422.04(2)EIAl=误差：1：0.1%2：0.8%。一般阶次越高，误差越大。6）求振型及画振型图 将各i代入有限元方程 ()200KM=，求 0i。02233Tiivv=现 令，有 21v=2223223124312 0 (1254)(613)0 (88)(6 13)(23)0(1254)(6 13)12 156(622)(6 13)(23)(622)(44)llllvllllll+=+取10.001841=222224312 0 (1254)(6 13)21.46 0 12.446.11 0 (88)(613)(23)(1254)(613)12 156(622)(613)(23)(622)(44)lllllllllll+=+2222 0 7.93