第三章 用有限元素法建立结构振动的数学模型

第三章 用有限元素法建立结构振动的数学模型

3．1 引言

【工程要求】：

对于简单的连续结构，如单件的杆、板、梁，可以建立结构振动的偏微分方程，但对于杆、板、梁组成的复杂结构，仍然采用建立偏微分方程的方法则十分困难。如果用假设模态法（李兹方法），对实际工程结构假设出品质良好的整个结构的假设模态也十分困难。要对结构振动进行数值分析，必须建立振动的数学模型——振动方程。

工程结构振动分析中，要采用将结构离散为有限自由度系统的方法——有限元素法，来建立结构的数学模型。

【发展简况】

有限元素法，是在上一世纪五十年代中期，经过M.T.Turner及J.H.Argyris等人的开拓性工作以及后来许多研究者的大量工作，发展起来的一种结构分析的有效方法，上一世纪六十年代初，由J.S.Archer及J.H.Argyris等人引入到结构动力学分析中来。

有限元素法发展到今天，已经非常成熟，而且与先进的计算机技术结合，已经形成了一个以有限元分析方法为基础的计算机辅助工程（CAE）的技术领域以及更进一步的虚拟产品设计（VPD）这样的先进概念。世界上著名的CAE分析软件商主要有MSC.software和Ansys等公司的产品。

【有限元动力学分析的任务】

在结构振动分析领域，有限元素法处理的问题主要是两类：结构固有振动特性计算和结构振动响应计算（包括频率响应分析与响应时间历程分析）。两类问题中，用有限元法建立振动数学模型是最基础的工作。

【有限元素法（分析结构振动问题）的特点】：

原则上，有限元素法由于其对复杂边界的适应性，它可以处理任何复杂的结构。

1

求解结果的精度可以根据需要不断改善，建模过程规范统一，计算形式适合于计算机求解。 【存在的问题】：

随着精度要求的不断提高，所要求的计算机容量和计算时间急剧增加，从而引出了大型特征值问题的快速求解方法、将大型结构振动问题转化为若干小型结构振动问题集合的子结构求解方法，以及结构振动问题的并行求解方法等问题的研究。

【工程结构振动分析方法】

从结构振动分析的发展历史看，经典的方法有：

1．集中质量法——将质量分别集中在若干节点处，形成集聚质量阵。结构的刚度仍然连续分布，采用材料力学中求柔度的方法，求出柔度系数，得到柔度矩阵，即用柔度法来形成刚度矩阵。

集中质量法存在问题：

对大型复杂结构，用材料力学的方法，进行柔度矩阵的求解显然是不现实的。

2．假设模态法——以李兹法为基础，选择一组假设模态组成的模态矩阵，对结构进行离散，如第二章所述的方法。

假设模态法存在的问题：

1．对几何形状复杂的结构，假设模态难以选择。

2．对整个系统用假设模态法得到的运动方程是高度耦合的，求解困难。 3．对不同的结构，要根据实际情况选取不同的假设模态，求解过程不规范统一。

引入有限元素法的思想既解决了上述方法的缺点，又保留了它们的优点。

【有限元法分析振动问题的基本原理】

2

用有限元法分析结构动态问题的基本思想，与结构静态分析的思想是一样的。它采用的方法仍然是：将结构分解为有限数目的单元，各元素间由节点相连，各单元内结构的变形用位移形函数（相当于元素级的假设模态）来表示，以节点位移作为控制变量（元素的广义坐标）。元素间的位移连续条件通过引入的形函数来满足，动态平衡条件通过最后导出的有限元方程来体现。由于节点数目是有限的，最后得到的方程是一个多自由度、离散的、线性的矩阵微分方程。

3．2 运动方程的建立

仍然采用熟悉的拉格朗日方程法建立其数学模型（运动方程）。 对任一单元内部任一点的位移{d}与节点位移{?}的关系：

{d}?[uvw]T?[N]{?}

e （3－1）

[N]称为假设的已知位移形函数（可以看成是单元的假设模态，一般仍采用静态变形函数） 显然：

?}?[N]{??e} {d （3－2）

单元的动能：

Te??1212?V?}{d?}dV??A{dT12eT{??}?V?A[N][N]dV{??}Te （3－3）

eTee{??}[m]{??}[m]?e?V?A[N][N]dV称为单元质量矩阵，质量阵是对称矩阵。

T

整个结构的动能为：

T??Te?e?e11eTee?}T[M]{q?} （3－4） {??}[m]{??}?{q22{q}是全结构的节点位移列阵，[M]??[mee]为全结构质量阵，?e代表对整个

结构各单元的组集（Assemble）。

单元的应变向量:

3

{?}?[B]{?}，{?}?[E]{?}

e （3－5）

[B]——几何矩阵，[E]——弹性矩阵，{?}——应力向量

单元的势能为：

Ue??121?2eVT{?}[E]{?}dV?eeT12{?}(?[B][E][B]dV){?}VeTTe（3－6）

{?}[k]{?}

全结构的势能：

U??Ue?e?2e1{?}[k]{?}?eTee12T{q}[K]{q} （3－7）

[K]??[kee] （3－8）

作用在单元上的分布力{f}的虚功：

?We??V{d}{f}dV?{?}TeT?V{N}{f}dV （3－9）

T单元节点力（广义力）

{pe}??V{N}{f}dV

T （3－10）

全结构的外力虚功：

?W???Weee?{q}{P}

T （3－11）

{P}??{pe} （3－12）

【阻尼的处理】

采用粘性阻尼假定：阻尼力与运动速度成正比，方向与速度相反。

单元中分布阻尼r的耗散函数（瑞利耗散函数）：

Re??1212?V?}{u?}dV?r{uT12eT{??}?VTer[N][N]dV{??} （3－13）

eTee{??}[c]{??}耗散力（即瑞利耗散力）与耗散函数的关系为：

{QD}??Ree?[c]{??e} （3－14） e?{??}i全结构的耗散函数：

4

R??eRe??e11eTee?}T[C]{q?} （3－15） {??}[c]{??}?{q22[C]??[cee] （3－16）

将全结构的动能、势能、耗散函数和广义力代入非保守系统的拉格朗日方程，

d?idt?q(?T)??T?qi??U?qi??D?i?q?Qi(i?1,2,?N)

（3－17）

得到：

??}?[C]{q?}?[K]{q}?{P} [M]{q （3－18）

【几个相关问题】：

1．进行实际结构的振动分析时，在各个单元的矩阵组集之前，还要对单元矩阵进行由单元的局部坐标系向结构的总体坐标系转换。

记局部坐标系下节点位移向量{?e}向总体坐标系下节点位移向量{?}的转换阵为[?e]，则坐标转换关系为：

{?}?[?][L]{?}

eee （3－19）

其中[Le]为对单元矩阵组集时“对号入座”的定位矩阵。

[m]?[L][?][m][?][L]

eTeTeee

（3－20） （3－21） （3－22）

[k]?[L][?][k][?][L]

eTeTeee[c]?[L][?][c][?][L]

eTeTeee[p]?[L][?][p]

eTeTe （3－23）

经过上述变换后的单元矩阵可以直接叠加得到结构总体矩阵。

2．为了求解结构的固有振动特性，需要求解无阻尼情况下结构的自由振动方程：

??}?[K]{q}?{0} [M]{q (3－24)

将固有振动的简谐运动形式

{q}?{?}sin?t

（3－25）

代入得到的结构的特征方程：

(??[M]?[K]){?}?{0}或 [K]{?}??[M]{?}

5

2 （3－26）

数学上构成所谓的广义特征值问题。

3．振动分析中采用的质量阵问题

在结构振动分析中，常采用的质量阵形成方法有：集中质量模型和一致质量模型。

在采用集中质量模型时，一般是按照杠杆原理将单元质量向单元各个节点上进行分配，在局部转动效应显著时，还要考虑单元的转动惯量。集中质量模型得到的质量阵为对角矩阵。

将单元内惯性分布视为与静力形函数同样规律的分布，导出的质量矩阵称为一致质量阵。即上面（３－２０）推导出的质量阵。这样的质量阵为满阵。

注意：集中质量阵与一致质量阵都不是振动结构在实际上精确的质量分布模型。理论上，结构的动位移是与频率相关的。动位移在不同振动频率和振型下是不同的，即不同的频率对应有不同的惯性。所以严格地讲，质量阵也是与频率相关的。下面以轴向振动的杆元为例，说明这个问题。

作为连续体的二力杆元的振动偏微分方程（波动方程）为：

x

E?d2??x2??d?t22?0 （3－27）

1 d1

2 L d2

记c?E/?称为波速，E，?分别为杆的弹性模量和密度。 方程（3－27）的解为：

d?[(cos?xc?cos?Lcsin?xc)csc?Lcsin?x?d1?]??ec?d2?i?t （3－28）

故此时形函数为:

[N]?[(cos?xc?cos?Lcsin?xc)csc?Lcsin?xc] （3－29）

?L代入单元质量阵公式得：

?L?L??Lcsc?cos?ALc?L?ccc[m]?csc??L?L2?Lc?1?cotcc???cc（3－30）

?L?L?L??csc?cosccc?1?cot?L显然，这个质量阵为满阵，且各元素为与频率相关的量。理论上这样的质量

6

阵真实反映实际的惯性分布，计算得到的固有特性更精确，但却使计算大大复杂化，而且只有对简单的情况（如二力杆元）下，才能从偏微分振动方程得出（动力）形函数，因此，工程实际中很少采用这种方法。而是采用与推导刚度阵时一致的（静力）形函数。

至于采用集中质量阵还是一致质量阵，得到的结果更好，没有固定的规律和结论，要视具体情况而定。从一般经验上讲，在单元划分较细时，用集中质量阵较好，反之宜采用一致质量阵。根据特征值隔离定理知道，采用一致质量阵分析得到的是固有频率精确解的上界，随着分元的细化，计算结果单调地向精确解逼近，而集中质量阵给出的结果就不具备这种特性，可能偏低也可能偏高。从工程分析经验看，由于建立有限元模型的离散过程，已经使结构比实际结构的刚度增大，因此采用集中质量阵可能有时反而会得到误差较小的结果，但这需要经验和技巧。但采用集中质量阵得到的振型一般误差较大，对振型要求较高时，还是宜采用一致质量阵。

集中质量阵是对角阵，在计算时可以节省计算时间。我们希望能获得一种优于集中质量阵的对角化质量阵。例如，对于梁的弯曲振动，可以按下式来计算对角化的非一致质量阵的单元：

m(i)rr?(??AN(y)dy)(??Ady)004li2rlili(r?1,2,3,4) （3－31）

??r?10?ANr(y)dy2而非对角元全部置零，Nr(y)(r?1,2,3,4)是梁的四个形函数。由于这种质量矩阵在一定程度上反映了单元形函数的特性，因此可以给出精度较好的计算结果。

3．3 典型结构单元的有限元建模

结构有限元模型的自由度数＝节点数×节点位移数

一、

纵向振动杆元

杆元i?j的任意截面处位移由两节点位移插值得到。

7

u(x,t)?Ni(x)ui(t)?Nj(x)uj(t)?[NiNiNNj]{uiuj}?[N]{q(t)}e （3－32）

j应满足边界条件：

Ni(0)?1,Nj(0)?0Ni(l)?0,Nj(l)?1 （3－33）

杆元受节点力时的静态方程为：

??x(AE?u?x)?0

（3－34）

xl其解为：

u(x)?c1?c2 （3－35）

xl代入边界条件得到形函数：

Ni(x)?1?xl,Nj(x)? （3－36）

从而杆元的质量阵和刚度阵为：

[m]???l0?A[N][N]dx?T?Al?2?6?11??2? （3－37）

?1??1?[k]?l0EA[N?][N?]dx?TAE?1?l??1 （3－38）

二、

弯曲振动梁元

梁元i?j的任意截面处位移由两节点位移插值得到。

y(x,t)?Nix(x)yi(t)?Ni?(x)?i(t)?Njx(x)yj(t)?Nj?(x)?j(t)?[NixNi?NjxNj??yi?????i?e]???[N]{q(t)}y?j???j???（3－39）

yi?iyj?j为梁元的节点位移。

均匀梁元受常值节点力作用时的挠曲线偏微分方程为：

?22?x(EI?y?x22)?0??y?x44?0 （3－40）

从而

8

y?c1?c2y(x)应满足边界条件：

xx2x3?c3()?c4() lll （3－41）

y(0)?yi,y?(0)??iy(l)?yj,y?(l)??j （3－42）

代入边界条件得到：

c1?yic2?l?ic3??3yi?3yj?2l?i?l?c4?2yi?2yj?l?i?l?jj （3－43）

代入位移表达式，整理后，得到形函数为：

?iy?1?3()?2()llxl2x2x3?i??x?2l()?l()lx3?jy?3()?2()llxl2x2x （3－44）

3?j???l()?l()lx3形函数阵：

[N]?[?iy?i??jy?j?] （3－45）

梁元的质量阵和刚度阵为：

?156??Al?22lT?A[N][N]dx?420?54???13l?12?6lEIT???EI[N][N]dx?3l??12??6l22l4l25413l1562[m]??l013l?3l6l4l2?22l?12?6l12?6l?13l?2??3l? （3－46） ?22l?2?4l?6l?2?2l? （4－47） ?6l?2?4l?[k]??l0?6l2l2上面给出的是平面梁元的特性矩阵，对于空间梁元特性矩阵形成过程完全相同。可以参考任何一本有限元素法方面的专著或教材。这里不再赘述。

【单元的坐标变换阵】：

单元局部坐标与总体坐标不一致时，需要进行坐标变换，记局部坐标X?Y

9

?ui???与总体坐标X?Y之间夹角为?，则局部坐标下节点位移列阵?vi?与总体坐标

????i??ui???下位移列阵?vi?间的变换关系为：

????i??ui??cos?????vi????sin?????0?i??sin?cos?00??ui??ui??????0?vi??[?]?vi? (3－48) ??????1????i??i?则坐标变换阵为：

?[?][?]???[0][0]?? [?]? (3－49)

三．面内振动的平板元

以三节点三角元为例，取三角元三个顶点的位移为单元广义坐标，

{q}?[uiviujvjumvm] （3－50）

eT单元内任一点的位移为：

?u(x,y,t)??Ni{e}??????v(x,y,t)??00NiNj00NjNm00?ee?q?[N(x,y)]{q(t)}（3－51） Nm???三个形函数Ni,Nj,Nm应满足边界条件：

?0Ni(xj,yj)???1i?j,i?j,i?i,j,mj?i,j,m （3－52）

设满足此条件的形函数为：

Ni(x,y)?ai1?ai2x?ai3yi?i,j,m （3－53）

代入边界条件可得到：

Ni(x,y)?12?(ai?bix?ciy)i?i,j,m （3－54）

??12111xixjxmyiyjymai?xjym?xmyibi?yj?ymci?xm?xj（3－55）

?为三角形单元的面积。aj,bj,cj,am,bm,cm由上式按下图轮换下标求得：

10

本文来源：https://www.bwwdw.com/article/h9k3.html