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

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

《动力学问题的有限元法.pdf》由会员分享，可在线阅读，更多相关《动力学问题的有限元法.pdf（22页珍藏版）》请在冰点文库上搜索。

第七章动力学问题的有限元法第七章动力学问题的有限元法结构动力学是研究动载荷作用下结构动力反应规律的学科，讨论结构在动力荷载作用下反应的分析方法，寻找结构固有动力特性、动力荷载和结构反应三者间的相互关系。

研究结构在动力荷载作用下的反应规律，能够为结构的动力可靠性（安全、舒适）设计提供依据。

前面介绍的静力学问题的研究对象是受不随时间变化的载荷作用。

而动力学问题的对象受随时间而变的载荷的作用，从而使在结构中产生的位移、速度、应力和应变都随时间而变。

当结构受随时间变化的载荷作用，且这种载荷的作用对结构的变形和应力的产生起主要作用，以致影响设备的安全性，或舒适性。

这时就要进行动力学分析，充分认识其规律性，从设计阶段就抑制这种不利状况的发生。

例如，有时虽然动载荷不大，但结构在交变力的作用下，其某些固有频率与激励力的作用频率相接近时，就会引起很大的振动、变形或应力，这时，就必须对结构作动力学分析。

又如，要利用结构在周期性作用力驱动下的定向振动，例如利用这种运动输送产品，这时，就必须巧妙地设计结构，使其具有某些与激励频率一致的固有频率，并且使结构对激励具有适当的响应能力。

总之，不管是利用振动，还是抑制振动，都需要进行结构动力学分析。

当前结构动力学的研究内容有三类。

第一类问题：

反应分析（结构动力计算），第二类问题：

参数（或称系统）识别，第三类问题：

荷载识别。

第一类问题是已知系统动态特性和动载荷作用部位及大小，求出系统的响应随时间变化的位移，速度，加速度和应力等。

第二类问题是已知系统的输入输出特性，分析系统固有的动态特性，结构模态分析就属于这一类问题。

第三类问题是在已知系统动态特性的条件下,通过测量系统的响应,或由响应准则预先给出响应要求,以此识别对响应的外载荷。

三类结构动力学研究内容的载荷、结构和响应之间的关系如图7-1所示。

动载荷种类大致分类如图7-2所示。

图7-1结构动力学研究的内容图7-2动载荷种类1本章主要介绍结构动力学分析的基础知识，并主要介绍系统固有特性的有限元分析方法有限元模态分析。

主要知识点：

1.导出有限元的结构动力学方程和模态分析有限元方程。

2.介绍模态分析方程中质量矩阵与阻尼矩阵的计算。

3.把计算结构的固有频率与振型的问题归纳为一个求特征值的问题。

7.1结构动力学方程及有限元方程7.1结构动力学方程及有限元方程系统的受迫振动的微分方程：

m（）xcxkxRt+=?

其中：

质量，阻尼系数，弹性系数。

mck,xxx?

分别为加速度，速度和位移。

1.单元的动力方程：

（用“虚功原理”求解）结构动力学有限元分析与静力学的有限元分析一样，首先要建立单元的有限元方程。

为求解起见，要引入达朗贝尔原理。

质点的达朗贝尔原理可表述为：

当非自由质点运动时，主动力，约束反力和惯性力构成一个动平衡力系。

惯性力是一种虚构的力，由质点的加速度引起。

mFNS即：

0fns+=现考察单元的动力学方程的建立。

1）.准备工作：

将位移函数表达成近似函数。

单元内任一点位移（,）（）eNNxyz=et，是单元内位移的近似函数,这儿N只是位置的函数，与时间无关。

而e与时间相关。

（,）（）eeNNxyzt=eeNN=?

2）.分析单元上的作用力。

图7-3单元分析主动力：

集中力，面力，体积力：

psRRRRV=+阻尼力：

与速度相关，并与速度方向相反，两者的比例关系为阻尼系数：

ceRCCN=?

2惯性力：

与加速度，质量相关：

单位体积惯性力meRN=?

阻尼力、惯性力与主动力方向相反。

密度，单位体积质量。

3）.建立动力学平衡方程：

利用虚功原理和达朗贝尔原理，系统处于动平衡，首先将载荷移置到节点上去。

利用虚功原理。

外力做的虚功：

11*11*TeTTppTeTTSsSsTeTTVVVVTeTTVCVCTeTTVmVmRNRRdSNRdSWRdVNRdVNRdVRdVNRdV=集中力：

面力：

体积力：

阻尼力：

惯性力：

RdV内力做的虚功：

*TeTTeVVUdVBDB=dV0虚功原理：

（）WUWU=，最小势能原理根据虚功原理，应有WU=1*1*eTeTeTTTpSsVVeTeTeTTTTeVCVmVNRNRdSNRdVNRdVNRdVBDBdV+=由*eT的任意性，可约去，且由阻尼力，惯性力的定义可知：

cmeeRCCNRN=?

，代入上式，有11TTTpSsVVeeTTTeVVVNRNRdSNRdVBDBdVCNNdVNNdV+=+?

=eRi令：

11TTTepSsVVTeVTeVTeVNRNRdSNRdVRKBDBdVcCNNdVmNNdV+=上式可表达为：

，这就是单元的动力学有限元方程。

eeemcK?

+=+=2.结构整体动力学有限元方程与静力学的集成原理相同：

根据局部单元节点编号与整体的节点编号的对应关系，利用，即一个节点上的节点力是该节点所在的所有单元的相应节点（同这个节点）上的节点力的总和来集成。

eeiiRR=3位移唯一性。

得MCK+=?

ReMm=，eCc=eKk=其中，当0C=，称为无阻尼受迫振动，此时MKR+=?

当0,0CR=，称为无阻尼自由振动，此时0MK+=?

7.2单元质量矩阵和阻尼矩阵的表达式7.2单元质量矩阵和阻尼矩阵的表达式在动力学有限元方程中除刚度矩阵之外，还有两个重要的系数矩阵：

即质量矩阵和阻尼矩阵。

这一节中将介绍单元质量矩阵和单元阻尼矩阵的表达式。

一.单元质量矩阵单元质量矩阵有二种表达形式：

一致质量矩阵和集中（堆积）质量矩阵。

1.一致质量矩阵：

用公式计算的单元质量矩阵，称为一致质量矩阵，是因其中的TeVmNN=dVN就是位移模式形函数N而得名。

或者说，推导质量矩阵采用的形函数与采用推导单元刚度矩阵时采用的形函数是一致的。

以平面应力问题为例：

三角单元，ijmNNININI=m（位移模式中的000000ijmijNNNNNNN=N）图7.3三节点三角形单元则三角形单元的一致质量矩阵em：

（以单位厚度计）ii22jj2j0N0N0N00N00N00jmijmijmeTijNNNNNNNNmNNdxdyNN=?

j2mm2mmNN0N000N0N0mijmijmdxdyNNNNNNN?

由于m，利用数值积分：

2iNii00N0N,iijjmNLNLNL=4!

2

（2）abcijmabcLLLdxdyabc=!

+?

，有22!

24!

62（112）!

12iijNdxdyNNdxdy=+从而得到20210212010210102010102em=对称2.集中质量矩阵：

将单元质量假想地集中到节点上去，也就是说只在节点上有质量。

这样，一个节点的加速度就不会影响其它节点的惯性力（相互之间不受影响）具体处理方法是将质量平均分到所有节点上去。

例如三角形三个节点单元，每个节点分到1/3的质量，即1m3it=,t为厚度，为三角形面积。

111111111331111ewmt=资料表明，集中质量矩阵计算出来的频率与一致质量矩阵的计算结果相差无几。

在单元相同的情况下，集中质量矩阵计算的结果低于一致质量矩阵的计算频率，而振型则是一致质量矩阵的计算结果精确。

由于集中质量矩阵式对角线矩阵，所以计算中常采用。

3.附加质量当有附加的质量在节点上时，要在计算整体质量矩阵时，加上这些附加质量cm即ecMmm=+二.单元阻尼矩阵根据阻尼的成因可以分为粘滞阻尼和结构阻尼。

1粘滞阻尼（比例阻尼）前面得到的阻尼矩阵是将阻尼看作是正比于节点的运动速度，称为粘滞阻尼。

它消耗振动的动能（能量）。

从可知，这个阻尼矩阵与质量矩阵成正比。

TeCCNNd=VdVeTmNN=eeCCm=，/Cc=如果,c均为常数的话。

所以称为比例阻尼。

（相对于质量矩阵而言）NTeCCNd=V2结构阻尼比例于结构的变形速度?

，由于材料内部的摩擦引起的，称为结构阻尼。

这时的阻尼力可以简5化为cRD=?

，这样可得到单元的结构阻尼矩阵使其具有的形式，所以结构阻尼正比于刚度矩阵：

BeTkCCBDd=VeekCk=3Rayleigh阻尼在实际分析中，精确决定阻尼矩阵是相当困难的，通常允许将整体阻尼矩阵简化为质量矩阵M和刚度矩阵K的线性组合，这样可用来解耦有阻尼的振动方程。

M+KC=，称为Rayleigh阻尼。

其中j22i222（-）2（-）ijiijjijjiji=其中22iii+=称为结构的第i阶振型的振型阻尼比振型阻尼比。

2i+称为结构的第i个校准振型的振型阻尼常数。

工程中一般取i=0.020.24，它与结构的类型，材料的能量消耗特性及振型有关。

可通过实测，从系统的固有频率和辐频特性曲线获得。

质疑质疑：

相关内容的矛盾描述：

在王勖成的书P476有如下字样：

“,是不依赖于频率的常数”。

但是，在同一本书的P487，又出现下列字样：

“如果根据实验或相近似结构的资料，已知二个振型的阻尼比i，j，可以得到j22i222（-）2（-）iiiijjijjiji=”7.3结构动力学特性的计算7.3结构动力学特性的计算求解结构的固有频率和对应于每一阶固有频率的振型，这是动力分析的基本内容，称为模态分析。

振型是指在某一个自振（固有）频率下，结构各节点振动值的集合，即节点振幅的相对大小。

由于阻尼对结构的自振频率和振型的影响不大，所以，一般只考虑无阻尼的自由振动。

一.模态分析的有限元方程在MPCK+=?

中，令P0,0C=，得模态分析的有限元方程：

6M0K+=?

（注：

注：

联系到温度场的瞬态热传导方程：

3PKtKt+=?

，P0=无内热源。

有人作比对分析，对温度场的分析采用所谓的“热模态分析”，但要注意，两个方程只是“形似”，它们的解有区别，热是一阶方程，振动是二阶方程（相对于时间）热是指数型解，振动是正弦函数分布。

）二。

数学预备知识二。

数学预备知识求解固有频率和振型，数学上涉及到矩阵的特征值问题。

这儿做一个有关矩阵的特征值，特征向量和特征方程方面的知识的扼要回顾，以便将我们的主要精力放到模态分析的物理意义上去。

1特征值和特征向量对n阶方程：

和n维非零列向量，如果有一个数1111nnnnnnaaAaa=?

12（,）Tn=?

使得矩阵nnA和维非零列向量n有A=，则称为矩阵A的特征值（特征根），向量为矩阵A的特征值所对应的特征向量。

（）0AAI=，为非零列向量，若等式成立，应有0AI=2.特征矩阵11121212211nnnnaaaaaaAIaan=?

称为特征矩阵，其中I为n阶单位矩阵3.特征多项式：

行列式AI称为矩阵A的特征多项式。

记为（）AI=4.特征方程式：

（）0=，称为A的特征方程。

特征方程（）0=的个根n12,n?

就是矩阵A的个特征值。

n集合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）成立，必须有002M0K+=即，要使0得到非零解，必须要有2M0K+=，从中可以求出i（i=1，n）,n为自由度数。

然后，将i代回到式

（2）求出每个i对应的0。

（一般工程上只需少数个最低的谐振频率即可），对应于每个ni，都有一组0i记2ii=广义特征值工程上称为谐振频率，自振频率。

0i特征向量工程上称为结构的振型，或振动模态。

2特征向量的性质1）0解不是唯一的，当0为解时，0也是解。

将0代入

（2）式时，有20（M）K+=00也是结构的一组特征向量，也就是说，在某个频率下，物体的各点的振动相互之间保持一个相对位移的比例关系，但它们的绝对值是可变的。

2）.0的归一化处理通常将0归一化处理，取0使得001TiiM=（3-1）这样归一化的0具有以下性质000010TiiTjiMM=（3-2）00000TiiiTjiKK=（3-3）证明：

由001TiiM=，有其余的表达式成立。

8证明有：

00TiiiK=证明过程:

由有限元方程0（）iiKM=0有00iiiKM=左乘得0Ti00000（）TTTiiiiiiiKMM=0i001TiiM=001=TiiiiK=证明有：

及000TjiM=当000TjiK=ji证明过程：

1）.对阶i0iiiKM0=a）对j阶00jjjKM=b）对a）式左乘0Tj有000TT0jiijKiM=

（1）对b）式左乘有0Ti000TTijjiK0jM=

（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方阵，nnX称为振型矩阵或模态矩阵。

9这时1001TnnXMXI=?

1200TXKX=?

，特征矩阵，频谱矩阵X是关于,MK的正交矩阵。

例一：

如图示悬臂梁，不计轴向位移，求其固有频率和振型。

图7-4悬臂梁1）位移模式无阻尼自由振动方程：

0MK+?

=,模态有限元方程：

（）200KM=特征方程为：

20KM=；不计纵向（轴向）位移，有：

ieijjjviv=，其中挠度,转角21234vaaxaxax=+3dvdx=边界条件：

110v=，所以223300evv=2）单元质量矩阵取两个单元。

单元质量矩阵eTVMNNdV=23232323232232132232xxxxxxxNxllllllll=+x102323232323232322322322323232132213223232elxxllxxxxxxxxxxxllMxAllllllllxxllxxll+=+dl+考虑到一维单元的自然坐标1ijxxLLLLL=，将eM转换成L表达的形式，可以使积分变得方便。

232232232232132（12）132（12）32（）32（）jjejjjjjjjjjjjjljjLLxLLMLLxLLLLxLLAdLLxLL+=+l+11112121111121211212222212122222xxyxxxyexxyyyxyyyxxyxxxyxxyyyxyyMMMMMMMMMMMMMMMMM=y23211（132）xxjjlMALLdl=+=1335lA（）23221111（132）（12）210yxjjjjlMALLxLLdllA=+=（）232321（132）（32）xxjjjjlMALLLLdl=+=970lA（）23221（132）（）yxjjjjlMALLxLLdl=+=213420lA（）2211（12）yyjjlMAxLLdl=+=21105lA23222（32）xxjjlMALLdl=1335lA（）2222（）yyjjlMAxLLdl=+=3105lA（）（）232121332（12）420yxjjjjl2MALLxLLdlAl=+=（）32212（）（12）140yyjjjjllMAxLLxLLdlA=+=（）（）2232211（）32210xyjjjjl2MAxLLLLdlAl=+=1122311112121222111121211212222222121222223322131191335210704201111321010542014091313117042035210131142014021010xxyxxxyxexyyyxyyyxxyxxxyxxyyyxyyyllllMMMMllllMMMMMAMMMMllllMMMMllll=52222156225413224133541315622420133224elllllllAMllllll=3）单元刚度矩阵梁的刚度矩阵可以用两种方法求得。

利用单元梁的平衡，利用能量法。

利用单元梁的平衡。

这儿要重温材料力学中关于梁的理论：

纯弯梁微分单元受弯矩M和剪力V的作用，梁的微分单元平衡时，有2323,=dvdvMEIVEIdxdx=。

图7-5梁的微分单元图7-6梁的有限元单元从图7-5和图7-6的比较中可以看出，,ijMMMM=；,ijVVVV=单元位移模式：

231234evaaxaxaxN=+=23232323232232132232xxxxxxxNxllllllll=+xieeeijjMVkRMV=，利用公式232,=dvdvMEIVEIdxdx=3，和位移模式evN=12可得33333232126126eeidNdvVEIEIEIdxdxllll=22222322321261266462eeiiiiidNxxxxdvMEIEIEIlldxdxllllll=+33333232126126eeidNdvVEIEIEIdxdxllll=22222322321261266462eejjjjjxxxxdNdvMEIEIEIlldxdxllllll=+现令0,ijxxl=，有323222223323222221261261261266462646212612612612662646264iiieeeijjjjjjllllvvVlliiMllllEIllllkREIvvVlllMllllllllllll=2232212612664621261266264ellllllEIKlllllll=，式中，l为单元长度，为单位长度质量。

利用能量法求单元刚度矩阵TeVKBDBdV=，几何矩阵B由几何方程求得。

eduBdx=022bdudxdvydx=是拉伸引起的正应变0dudx=和由弯曲引起的正应变22bydydx=v的总和。

本题不计轴向位移。

00dudx=2222ebdNdvyydxdx=22dNBydx=由evN=和23232323232232132232xxxxxxxNxllllllll=+x有22232236124661226dN2xxxxxByylldxllllll=+13232223223223261246612466122661226TeVVxllxxxxxllKBDBdVyEdVxllllllllllxll+=+22323223232322223222322612612466126126122646612464661246eVxxxxxxxllllllllllllllxxxxxxllllllllllllKEy+=2223232322323222322222661261261246612612262661226462661xllxxxxxxxllllllllllllllxxxxxlllllllllll+232226dVxxlll+22323223232322223222322612612466126126122646612464661246exxxxxxxllllllllllllllxxxxxxllllllllllllKEydA+=2223232322323222322222661261261246612612262661226462661xllxxxxxxxllllllllllllllxxxxxlllllllllll+232226Dldlxxlll+其中抗弯刚度。

2DEydAEI=令2232322323232222322232612612466126126122646612464661246xxxxxxxllllllllllllllxxxxxxllllllllllllF+=222323232232322232222326612612612466126122626612264626612xllxxxxxxxllllllllllllllxxxxxxllllllllllll+2226ldlxll+对其中的每一个元素分别积分，得到：

2112361212lxFd3lll=+=l，122323612466lxxlFdlllll=+=l，132323361261212lxxFdllllll=+=14232612266lxxFd2lllll=+=l，2222464lxFdllll=+=，232233466126lxxlFdllllll=+=14142246262lxxFdllllll=+=，23323361212lxFdlll=l，342322612266lxxFdllllll=+=2442264lxFlll=+=dl其他各元素可利用对称性求得。

2232212612664621261266264ellllllEIKlllllll=4）质量矩阵和刚度矩阵的集成Mmm=+，Kkk=+222223221261260064620012624012662086200126126006264llllllllEIKllllllllllll=，22222221562254130022413300541331205413420133081330054131562200133224lllllllllAMlllllllllll=（注：

注：

当有较多位移边界条件已知时，可以在集成计算各元素之前，先作降阶处理，划行划列，划去相应的行与列，只集成剩余的各元素。

这样可节约大量时间。

）5）求固有频率代入边界条件：

110v=，所以，22223120541308133541315622420133224lllllAMllllll=，2232224012608621261266264llllEIKlllllll=特征方程20KM=22222322222401263120541308620813301261265413156224206264133224llllllllEIlAlllllllllllll=令：

42420lAEI=则方程化为2222243120（1254）（613）0（88）（613）（23）0（1254）（613）12156（622）（613）（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=,2222.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=22232231243120（1254）（613）0（88）（613）（23）0（1254）（613）12156（622）（613）（23）（622）（44）llllvllllll+=+取10.001841=2222243120（1254）（613）21.46012.446.110（88）（613）（23）（1254）（613）12156（622）（613）（23）（622）（44）lllllllllll+=+222207.93