有限元第六章--动力问题的有限元法

1第六章动力问题的有限元法6.1概述前面几章所研究的问题都属于静力问题，其特点是施加到结构上的外载荷不会使结构产生加速度，且外载荷的大小和方向不随时间变化，因而结构所产生的位移和应力也不随时间变化。本章将要研究结构分析中另一类重要问题的有限元解法，即动力问题的有限元解法。动力学问题的特点是，载荷是随时间变化的，因而结构所产生的位移和应力是时间的函数，结构会产生速度和加速度。由于结构本身的弹性和惯性，结构在动力载荷的作用下，往往呈现出振动的运动形态。结构振动是工程中一个很普遍很重要的问题。有些振动对我们有利，例如，振动打桩，振动选料，有些振动对我们有害，例如，机床的振动，仪器与仪表的振动，桥梁、水坝及高层建筑在地震作用下的振动等。因此，我们必须对振动体本身的振动特性以及它对外部激振力的响应有一个明确的认识，才能更好地利用它有利的一面，而避免它有害的一面，设计出更好的机械和结构。振动问题主要解决两方面的问题。1.寻求结构的固有频率和主振型，从而了解结构的固有振动特性，以便更好地利用或减少振动。2.分析结构的动力响应特性，以计算结构振动时动应力和动位移的大小及其变化规律。6.2结构的振动方程结构的振动方程可用多种方法建立，这里我们使用达朗伯原理（动静法），仿照前几章建立静力有限元方程的方法，来建立动力问题的有限元方程。在静力问题中用有限元法建立的平衡方程是}{}]{[FK在振动问题中，对结构的各节点应用达郎伯原理所建立的振动方程仍然具有与上式相同的形式，只不过节点位移是动位移，节点载荷是动载荷，它们都是时间的函数。上面的方程成为)}({)}(]{[tQtK(6.1)上式中)(t为节点的动位移，它是时间的函数，)}(]{[tK是t时刻的节点位移产生的弹性恢复力，它与该时刻的节点外力)(tQ构成动态平衡。在动态情况下，结构承受的载荷（集中载荷，分布载荷）可随时间而变化，是时间的函数。按有限元方法将此种载荷移置到节点上，得到的节点载荷向量)}({tF也是时间的函数。此外，结构在运动中，各点除位移f以外，还有速度.f及加速度..f。按照达郎佰原理，有加速度的质量应附加有惯性力载荷。如材料的密度为，则结构单位体积的惯性力为}{..f。这对结构来说，相当于又受有另一种体积力，大小与点的加速度成比例，而方向与加速度方向相反。另外，在结构运动过程中，还会受到周围介质和来自内部的阻力。精确地描述这种阻力的变化规律是很困难的，一般采用阻力与速度f成比例的近似线性假定，如阻力系数为μ，则单位体积的阻力为}{f。这对结构来说相当于另一种体积力，大小与点的速度成比例，方向与速度方向相反。按有限元方法，用单元节点位移e进行插值表示单元内部位移。eNf}]{[(6.2)2此处形函数仍只是位置的插值函数，与时间无关，则单元内的速度和加速度分别为eNf}]{[(6.3)以及eNf}]{[(6.4)其中e、e为单元节点的速度及加速度向量。将单元惯性力f与阻力f作为体积力，按式(3.28)移置到单元各节点，就得到相应的单元等效节点载荷向量，记为eF和eF，则有dvfNFTe}{][dvfNFTe}{][将式（6.3）、（6.4）代入上式，有eeTeTemdvNNdvNNF][][][][][eeTeTecdvNNdvNNF][][][][][式中dvNNmTe][][][(6.5)称为单元质量矩阵。由于推导式（6.5）时采用了与推导单元刚度矩阵时相一致的形函数，故式（6.5）所表示的质量矩阵也称为一致质量矩阵。而dvNNCTe][][][(6.6)称为单元阻尼矩阵，由式（6.5）和（6.6）可见单元质量矩阵和单元阻尼矩阵是对称的。将移置到节点上的动载荷、惯性力、阻力作为载荷，按单元叠加，得到有限元节点位移方程：meeeeecmFK1)}{][}{]([}{}]{[或}{}]{[}]{[}]{[FKCM(6.7)其中meemM1][][(6.8)称为结构质量矩阵或总质量矩阵，而meecC1][][(6.9)称结构阻尼矩阵。可见结构质量矩阵和结构阻尼矩阵分别为单元质量矩阵和单元阻尼矩阵的叠加，其叠加方法与结构刚度矩阵的形成完全一样，借助于单元定位向量，用单元集成法完成叠加过程。由于em][和ec][是对称的，因而叠加合成的M和C也是对称的。3式（6.7）是节点位移的二阶微分方程，称为结构的动力方程式。对于不同的结构，可以选用不同的单元，有不同的形函数矩阵N，但动力方程（6.7）的建立过程都是一样的。当结构不受外载荷时，0F，如果再忽略阻尼，则动力方程（6.7）式成为0}]{[}]{[KM(6.10)这是系统的自由振动方程。弹性结构的振动实际上是连续体的振动，位移f是连续的，具有无限多个自由度。经有限元离散化处理后，单元内的位移按假定的位移形式来变化，可用节点位移插值表示。这样，连续系统的运动就离散化为有限个自由度系统的运动了。如果全部节点有N个自由度，则式（6.10）就是N阶的自由振动微分方程了。6.3一致质量矩阵与集中质量矩阵按照式（6.5），当单元的位移插值形函数矩阵确定后，就可用此式算出它的一致质量矩阵。例如，对平面三节点三角形单元，按照2.3节的叙述，形函数矩阵N可用面积坐标表示为][][mjiILILILN其中I为二阶单位阵。将上式代入式(6.5)可得平面三节点三角形单元的一致质量矩阵dxdytLILLILLILLILLILLILLILLILLILdxdytILILILILILILdxdytNNmmmjmimmjjjijmijiiimjimjiTe][][][][利用积分公式(2.56)，可由上式求得2104104100210410414102104100410210414104102100410410213][Mme(6.11)用式(6.5)可以计算出其它类型单元的质量矩阵。例如平面桁架单元的质量矩阵为421616131][Mme(6.12)而平面刚架单元的质量矩阵为LLLLLLLMme42203130156013540140007042201560140420][22称对(6.13)式(6.11)、(6.12)及(6.13)中的Ｍ为单元质量。由单元一致质量矩阵叠加形成的结构质量矩阵，一般都是稀疏、带状的，但都有相当的半带宽，如果将单元质量矩阵近似作为对角型矩阵，单元叠加后的结构质量矩阵也是对角阵。这种对角型的质量矩阵称为集中质量矩阵，而对角型的集中质量矩阵对结构的动力分析是非常有利的。可以将单元的质量以某种方式分配在单元的节点上而得到单元集中质量矩阵。例如当质量均匀分布时，将质量平均分配给各节点。按照这种方法平面三节点三角形单元的集中质量矩阵为][3110110113][IMMme(6.14)其中I为６阶单位阵。Ｍ为单元质量。而四节点四面体单元的集中质量矩阵为][4][IMme(6.15)其中[I]为12阶的单位阵。对于矩形弯曲板单元，在不计转动质量影响的条件下，它的集中质量矩阵为5001000100010014][Mme(6.16)对于平面刚架单元，在不计转动质量影响的条件下，它的集中质量矩阵为010100112][Mme(6.17)6.4阻尼矩阵各种工程结构的阻尼力及其产生的机理是非常复杂的。从宏观上看，阻尼有两种主要形态。一种是结构周围粘性介质产生的阻尼，称为粘性阻尼。粘性阻尼的阻尼力一般近似与运动速度成正比。另一种是结构材料内部磨擦产生的阻尼，称为结构阻尼或材料阻尼。结构阻尼的阻尼应力一般近似地认为与弹性体的应变速率成正比。如果假定阻尼力与运动速度成正比，那么在运动的弹性体中任意点处单位体积上作用的阻尼力为][}{NftPd式中——比例常数；——材料密度；N——单元形函数矩阵；t——单元节点速度矢量。可以将阻尼力}{dP看成是一种体积力，其等效的单元节点阻尼力向量为6vTvdTeddvNNdvPNF}{][][}{][}{或写成}]{[}{cFed因此，单元阻尼矩阵veTemdvNNC][][][][(6.18)它正比于单元质量矩阵em][。如果假定阻尼应力与弹性体的应变速率成正比，则阻尼应力可表示为edBDtD}]{][[}{][式中，为比例系数，D为弹性矩阵，B为应变矩阵，e为单元节点速度向量。下面推导阻尼应力的单元等效节点阻尼力向量。设单元等效节点阻尼力向量仍用edF表示，设节点发生虚位移}{*，单元内各点产生的虚应变为}{*，则}]{[**B。edF}{在虚位移}{*上做的虚功为edTeFW}{}{*单元的虚应变能为vvdTTdTedvBdvU}{][}{}{}{**由eeUW得到evvTdTeddvBDBdvBF}{]][[][}{][}{或写成eeedcF}{][}{因此，单元阻尼矩阵veTekdvBDBc][]][[][][(6.19)它正比于单元刚度矩阵ek][在实践中，要精确地计算各种单元的阻尼矩阵是很困难的。通常在程序设计中，假定结构总体阻尼矩阵C是结构总体刚度矩阵K与总体质量矩阵M的线性组合。称为瑞利阻尼，其表达式为][][][KMC(6.20)其中，比例系数及可通过实验确定。7采用瑞利阻尼近似，可以使运动方程求解大大简化，并且在程序中不必单独存贮总体阻尼矩阵。在实际工程问题中，阻尼的作用对结构的动力响应的影响并不大，这种近似处理具有实用价值。6.5结构的自振特性和特征值问题结构的自振特性是指结构的振动频率和振型，求结构的自振频率和振型也称对结构进行模态分析，是结构动力计算的主要内容之一。计算经验指出，结构的阻尼对结构的频率和振型的影响很小，所以求频率振型时可以不考虑阻尼的影响。此时系统的自由振动方程如式(6.10)，即0}]{[}]{[MK(6.21)当系统作自由振动时，各质点作简谐振动，各节点的位移可表示为tcos}{}{(6.22)将(6.22)代入(6.21)，并消去公因子tcos得到0}]){[]([2MK或}]{[}]{[2MK(6.23)因此，求解(6.21)式就是寻找满足式(6.23)的2值和非零向量，这种问题称为广义特征值问题。记2，和分别称为广义特性值和广义特征向量。式(6.23)可以写成0}]){[]([MK(6.24)这是一个齐次的线性方程组，若要有}{的非零解，系数行列式必须等于零，即0][][MK(6.25)展开此式可得0221122222221221112121111nnnnnnnnnnnnMKMKMKMKMKMKMKMKMK如果弹性结构的总体刚度矩阵K和总体质量矩阵M的阶数是n，则上述行列式展开后成为的n次代数方程式，由此可以求出n个根，即n个广义特征值nii,..