有限单元法(第十章)

结构物可能受到强烈地震、机械振动或海浪冲击等动荷载作用，使结构发生不允许的变形甚至破坏。因此，结构设计必须考虑动荷载效应，进而要求进行动力学分析。河北工业大学土木工程学院第十章动力分析本章以动力学问题为例，阐明如何应用有限元进行结构动力分析，内容包括：（1）动力有限元方程及系数矩阵；（2）结构固有特性计算；（3）结构动力响应计算；（4）解的稳定性分析。河北工业大学土木工程学院10.1动力有限元方程10.1.1有限元方程第十章动力分析振动系统中由质量、阻尼和刚度所引起的三种力以及外荷载，将振动系统离散后，便可得到有限元动力平衡方程组。（1）一般方程（10.1）其中，M、C、K分别为整体质量矩阵，阻尼矩阵；分别为整体节点加速度，速度和动荷载向量。河北工业大学土木工程学院10.1动力有限元方程10.1.1有限元方程第十章动力分析对于地震反应问题，直接作用在结构上的荷载P=0,且阻尼力和弹性恢复力分别只与相对速度和相对位移有关，而惯性力则与绝对加速度有关。现仍以分别表示相对加速度、相对速度和相对位移，表示地震牵连加速度（地面加速度）则式（10.1）变为（2）地震问题（10.2）其中（10.3）河北工业大学土木工程学院10.1动力有限元方程10.1.2方程的推导第十章动力分析在动力分析中，同静力分析一样首先对空间进行离散。单元位移函数为（1）位移、速度、加速度（a）为简便，我们把动力问题化成静力问题，即把惯性力和阻尼力看做体积力施加在结构物上，这样就可按静力问题进行分析。或河北工业大学土木工程学院10.1动力有限元方程第十章动力分析单元内部各点的速度和加速度分别可以用单元节点的速度和加速度表示为其中节点位移向量为（a）10.1.2方程的推导形函数矩阵为（b）河北工业大学土木工程学院10.1动力有限元方程第十章动力分析由式（1.19）惯性力的等效节点荷载（2）惯性力的等效节点荷载（d）若材料的质量密度为ρ，根据D’Alembert原理，单元内单位体积上作用的惯性力为其中，Me为单元质量矩阵，即（10.4）（10.5）10.1.2方程的推导河北工业大学土木工程学院10.1动力有限元方程第十章动力分析假设阻尼力正比于运动速度，比例常数为α，则单元内单位体积上作用的阻尼力可表示为（3）阻尼力的等效节点荷载（f）阻尼力的等效节点荷载为（g）其中，Ce阻尼力的等效节点荷载为10.1.2方程的推导（10.7）河北工业大学土木工程学院10.1动力有限元方程第十章动力分析用P表示外部动荷载产生的整体节点荷载向量，则总的整体节点荷载向量为。（4）动力有限元方程的建立（10.8）动力有限元方程式为10.1.2方程的推导河北工业大学土木工程学院10.1动力有限元方程10.1.3系数矩阵第十章动力分析采用有限单元法进行结构动力计算，必须建立结构系统的整体质量矩阵M、阻尼矩阵C和刚度矩阵K。（1）质量矩阵（10.1）整体质量矩阵M的元素mij称为质量影响系数，其物理意义是：自由度j的单位加速度在自由度i方向引起的力。在对结构进行离散化处理时，分配单元质量的常用方法有两种，即一致质量法和集中质量法。河北工业大学土木工程学院10.1动力有限元方程第十章动力分析对于平面问题常应变三角形单元，一致质量矩阵为一致质量法按式（10.5）计算单元质量矩阵。此时使用了与推导单元刚度矩阵相同的位移函数，故称为一致质量矩阵。又因为单元的动能和势能是相互协调的，也称为协调质量矩阵。10.1.3系数矩阵河北工业大学土木工程学院10.1动力有限元方程第十章动力分析关于质量矩阵的对角化，最常用的方法是对一致质量矩阵的行求和，即集中质量法简单地将单元的质量集中分配于单元的节点，形成集中质量矩阵。假定质量集中在点上，一般不考虑转动惯量，所以与转动自由度有关的质量系数为零。集中质量矩阵是对角矩阵。10.1.3系数矩阵

河北工业大学土木工程学院研究表明，采用一致质量矩阵将高估自振频率；而采用集中质量矩阵则以同样量级低估最高自振频率。因此有学者建议在实际中采用混合质量矩阵。即这两种矩阵的平均值。10.1动力有限元方程10.1.3系数矩阵对于平面问题常应变单元，集中质量矩阵为（10.10）第十章动力分析

河北工业大学土木工程学院动力学方程的阻尼项代表在运动中所消耗的能量。产生阻尼的原因是多方面的，如滑动摩擦、空气阻力、材料内摩擦等。通常是用等效粘滞阻尼来代替。10.1动力有限元方程10.1.3系数矩阵（2）阻尼矩阵（10.7）假设阻尼力正比于运动速度，则单元阻尼矩阵正比于单元质量矩阵，如式（10.7）。第十章动力分析

河北工业大学土木工程学院于是，根据初应力的等效节点荷载公式（9.50），可以得到单元阻尼矩阵10.1动力有限元方程10.1.3系数矩阵如果假定阻尼应力正比于应变速度（由于材料内摩擦引起的结构阻尼通常可以简化为这种情况），则可表示为（9.50）（10.11）第十章动力分析

河北工业大学土木工程学院10.1动力有限元方程10.1.3系数矩阵可见，此单元阻尼矩阵正比于单元刚度矩阵。

Rayleigh阻尼为质量矩阵和刚度矩阵的线性组合。（10.12）第十章动力分析（10.11）

河北工业大学土木工程学院10.2.1广义特征值问题10.2结构固有特性计算结构的固有频率和固有振型称为特征值问题。它是结构动力分析的基本内容，也是采用振型叠加法计算结构动力响应的前提。（10.1）在式（10.1）中令P=0，得到自由振动方程。在实际工程中，阻尼对结构自振频率和振型的影响不大，因此可忽略阻尼力，从而得到无阻尼自由振动的运动方程（10.13）第十章动力分析

河北工业大学土木工程学院将其代入式（10.13），可得齐次方程在自由振动时，因结构中各节点的振幅不全为零，故其系数行列式必为零，即（10.14）式（10.13）为常系数线性齐次常微分方程组，其解的形式为10.2.1广义特征值问题10.2结构固有特性（10.15）（10.16）求解方程组（10.15）的问题称为特征值问题。第十章动力分析

河北工业大学土木工程学院对应每个特征值ωi2，由式（10.15）可确定一组相应的振幅值φi2(i=1,2,┄n），称为特征向量，在工程上通常称为结构的振型。10.2.1广义特征值问题10.2结构固有特性（10.16）K和M都是n阶方阵，n是结构节点自由度数目。可见，式（10.16）是ω2关于的n次代数方程，称为式（10.13）的特征方程。其解ωi2(i=1,2,┄n）称为特征值。第十章动力分析

河北工业大学土木工程学院数学上可以证明，当质量矩阵为对称正定、刚度矩阵为对称正定或半正定时，所有特征值为非负的实数，特征向量也是实向量。10.2.2特征值和特征向量10.2结构固有特性（10.17）在结构分析中，刚度矩阵和质量矩阵是实对称矩阵。消除刚体位移后的刚度矩阵是正定的；采用一致质量矩阵时，质量矩阵是正定的。其中，ωi称为结构的第i阶固有频率，对应的特征向量φi称为第i阶固有振型。第十章动力分析

河北工业大学土木工程学院显然，如果是广义特征问题的特征向量，乘以不等于零的常数后仍为特征向量。通常使特征向量满足10.2.2特征值和特征向量10.2结构固有特性（10.18）它们满足方程（10.19）这样规定的特征向量或振型称为正则振型。将特征解（ωi2,φi）和（ωj2,φj）分别代回方程（10.15）得到（10.20）（10.21）第十章动力分析

河北工业大学土木工程学院可见，当ωi≠ωj时，必有10.2.2特征值和特征向量10.2结构固有特性（10.22）式（10.20）两端前乘以φjT

，式（10.21）两端前乘以φiT

，并注意到K和M的对称性，可以得到上式表明特征向量关于矩阵M是正交的。由式（10.19）和(10.22)得（10.23）第十章动力分析

河北工业大学土木工程学院如果定义特征矩阵10.2.2特征值和特征向量10.2结构固有特性（10.24）由式（10.18）和（10.23）得（10.25）（10.26）（10.27）则特征解的性质可表示为第十章动力分析

河北工业大学土木工程学院振型叠加法的基本思想是：先将方程组非耦合化，然后再积分求解非耦合方程组。10.3.1振型叠加法10.3结构动力响应动力有限元方程组的解法主要有两类，即振型叠加法和直接积分法。M和K是振型正交的，而且在采用系统比例阻尼假设下，C也是振型正交的。所以借助振型向量所组成的位移变换矩阵，便可实现方程组的非耦合化。第十章动力分析

河北工业大学土木工程学院其中，xi(t)是时间的函数，这样便有（10.28）在求出结构无阻尼自由振动的频率和振型以后，可用振型的线性组合来表示结构的节点位移10.3结构动力响应10.3.1振型叠加法（10.29）代入方程（10.1），并各项前乘ΦT，得（10.1）（10.30）第十章动力分析

河北工业大学土木工程学院于是（10.31）注意到式(10.26),并设阻尼矩阵为质量矩阵与刚度矩阵的线性组合，即10.3结构动力响应10.3.1振型叠加法（10.32）其中（10.33）第十章动力分析（10.26）

河北工业大学土木工程学院将式（10.26）,（10.32）代入式（10.30），可得一组相互独立的微分方程（10.34）显然，如果已知结构的两个频率ωi，ωj以及相应的阻尼比λi,λj，则阻尼常数为10.3结构动力响应10.3.1振型叠加法第十章动力分析

河北工业大学土木工程学院Duhamel积分基本思想是：将激振力ri(t)分解为一系列微冲量的连续作用，分别求出每个微冲量的响应，然后将所有微冲量的响应叠加起来。由于高阶振型对结构动力响应的贡献一般都很小，通常只计算最低的3阶到5阶振型即可。（10.35）10.3结构动力响应10.3.1振型叠加法振型叠加法的缺点是必须求解特征值问题。另外，由于该法应用了叠加原理，故只适用于线性问题。其中，ai,bi由初始条件确定。第十章动力分析

河北工业大学土木工程学院（2）在一定数目的Δt区域内，假设位移at，速度,加速度的函数形式,计算t+Δt时刻运动状态的公式直接积分法是：直接对方程（10.1）进行逐步积分，而不进行任何形式的变换。10.3结构动力响应10.3.2直接积分法第十章动力分析求解的基本思路基于如下两个概念：（1）将在0<t<T内任何时刻都应满足运动方程，代之以仅在相隔Δt的离散时间点上满足运动方程；对t的离散方法不同，就有不同的数值积分方法，如中心差分法、线性加速度法、Newmark法和Wilson法。

河北工业大学土木工程学院(1)中心差分法在以下讨论中，假定求解域被等分为n个时段，且在0,Δt,2Δt┄时刻的解已经求得，计算的任务是求t+Δt时刻的解。（10.36）10.3结构动力响应10.3.2直接积分法第十章动力分析在中心差分法中，加速度和速度可以用位移表示为为了求得t+Δt时刻的位移解答at+Δt，可在t时刻建立运动方程（a）（b）

河北工业大学土木工程学院其中，称为有效刚度矩阵；称为有效节点荷载。将式（a）和（b）代入上式，得到（10.37）10.3.2直接积分法第十章动力分析10.3结构动力响应（10.38）（10.39）（10.40）当t=0时，为了计算aΔt，由（a）和（b）式得到河北工业大学土木工程学院这种方法的优点是计算简单，缺点则在于它是有条件稳定的，即当时间步长Δt过大时，积分是不稳定的。10.3.2直接积分法第十章动力分析10.3结构动力响应因at+Δt是利用t时刻运动方程得到的，K不出现在刚度矩阵中，因此中心差分法称为显示积分法。式（10.37）是用相邻时刻的位移表示的代数方程组，由此可解出at+Δt。

10.3.2直接积分法第十章动力分析10.3结构动力响应

10.3.2直接积分法第十章动力分析10.3结构动力响应

10.3.2直接积分法第十章动力分析10.3结构动力响应

10.3.2直接积分法第十章动力分析10.3结构动力响应

10.3.2直接积分法第十章动力分析10.3结构动力响应

河北工业大学土木工程学院于是，t+τ时刻的速度和位移可积分求得（2）线性加速度法（a）10.3.2直接积分法第十章动力分析10.3结构动力响应假定在时刻[t,t+Δt]内加速度呈线性变化，即（b）（c）将式（a）代入上式，τ=Δt得到

河北工业大学土木工程学院由式（c）得到（d）10.3.2直接积分法第十章动力分析10.3结构动力响应将式（d）代入式（b）得到（e）t+Δt时刻的位移解答at+Δt

可通过满足t+Δt时刻动力学方程（10.41）将式（d）（e）代入式（10.41）到由

河北工业大学土木工程学院其中（10.42）10.3.2直接积分法第十章动力分析10.3结构动力响应将式（d）和（e）代入式（10.41））得到由计算的公式（10.43）（10.44）在线性加速度法中，at+Δt是利用t+Δt时刻的时刻运动方程得到的，K出现在刚度矩阵中。因此称为隐式积分法。可以证明，线性加速度法也是有条件稳定的。

河北工业大学土木工程学院其中，γ和β是按积分稳定性要求而决定的参数。当γ=1/2,β=1/6时，式(a)和(b)相应于线性加速度法。（3）Newmark法（a）10.3.2直接积分法第十章动力分析10.3结构动力响应Newmark（1959）假设（b）在Newmark法中，

t+Δt时刻的位移解答at+Δt

也是通过满足t+Δt时刻动力学方程得到（10.45）

河北工业大学土木工程学院代入式（a），再一并代入式（10.45），从而得由计算的公式，首先从式（b）中解得（10.46）10.3.2直接积分法第十章动力分析10.3结构动力响应其中（10.47）（10.48）

河北工业大学土木工程学院（4）Willson法10.3.2直接积分法第十章动力分析10.3结构动力响应积分上式可得（a）Willson法假定加速度在[t,t+θΔt]内呈线性变化，参数θ≥1.0。以τ表示时间增量，0≤τ≤θΔt，则在[t,t+θΔt]内有（b）（c）

河北工业大学土木工程学院在式（b）和（c）中令τ=θΔt，得10.3.2直接积分法第十章动力分析10.3结构动力响应由此可得到用at+θΔt表示的（d）（f）（e）（g）

河北工业大学土木工程学院将式（f）和（g）代入上式，得10.3.2直接积分法第十章动力分析10.3结构动力响应其中（10.49）（10.50）t+θΔt时刻的动力学方程为（10.51）（10.52）

河北工业大学土木工程学院由此式计算t+θΔt时刻的加速度、速度和位移。由式（10.50）解出位移at+θΔt，由式（f）确定。令τ=Δt，可得10.3.2直接积分法第十章动力分析10.3结构动力响应（10.53）河北工业大学土木工程学院在选择时间步长Δt时，必须考虑解的稳定性和精度。10.4.1基本概念第十章动力分析10.4解的稳定性无条件稳定：在任何时间步长Δt下，对于任何初始条件，方程的解不无限止地增长，则称此算法是无条件稳定。有条件稳定：如果Δt必须小于某个临界值Δtcr时，上述性质才能保持，则此算法是有条件稳定。河北工业大学土木工程学院由于各振型的运动是相互独立的，方程是相似的，故只须从中取出典型的运动方程来分析即可10.4.1基本概念第十章动力分析10.4解的稳定性解的稳定性实质上是误差的响应问题，因此可令ri=0。阻尼对解的稳定性是有利的，令阻尼项为零，可用下式讨论解的稳定性（10.54）（10.55）河北工业大学土木工程学院（1）中心差分法10.4.2稳定性条件第十章动力分析10.4解的稳定性对式（10.55）进行积分，根据式（10.37），可以写出（10.56.3）（10.56.4）假定解的形式为将上式代入（10.56.3),则可的特征方程（10.56.5）（10.56.6）（10.55）

河北工业大学土木工程学院λ的根关系到解的性质，首先为使在小阻尼情况下的解具有振荡特性，λ必须是复数，这就要求解出上式的根第十章动力分析因为，所以从上式得到10.4解的稳定性10.4.2稳定性条件（10.56.7）（10.56.8）（10.56.9）

河北工业大学土木工程学院式（10.56.7）表示λ1,2的，表示无阻尼的自由振动。其次为使解不会无限地增长，还应要求第十章动力分析直接积分法相当于利用同样的时间步长对所有n个振型的单自由度方程同时进行积分，因此中心差分法为保持解的稳定性，时间步长必须满足（10.56.11）10.4解的稳定性10.4.2稳定性条件（10.56.10）其中Δtcr是临界时间步长，Tn是系统的最小固有周期。

河北工业大学土木工程学院将Newmark方法的循环计算公式（10.46）用于式（10.55）表示的运动方程，可以得到（2）Newmark方法由下面两式可得（10.57.13）第十章动力分析10.4解的稳定性10.4.2稳定性条件（10.57.12）（10.57.13）

河北工业大学土木工程学院利用上式和（10.55），式（10.57.12）可以写成其中第十章动力分析10.4解的稳定性10.4.2稳定性条件（10.57.14）代入上式可以得到关于λ的特征方程仍假定解的形式为（10.57.15）

河北工业大学土木工程学院其中该方程的根是现在来分析界的稳定条件：第十章动力分析10.4解的稳定性10.4.2稳定性条件（10.57.16）1.真正的解在小阻尼情况下必须具有振荡的性质，因此λ应是复数，这就要求

河北工业大学土木工程学院当Pi很大时，即Δt不受限制时，仍要求上式成立，必须是亦即（2）稳定的解必须不是无限增长的，因此第十章动力分析10.4解的稳定性10.4.2稳定性条件（10.57.18）（10.57.19）亦即（10.57.20）

河北工业大学土木工程学院因为当条件（10.57.19）满足时，式（10.57.22）恒成立。所以Newmark方法无条件稳定的条件是同样，当pi很大时，仍要求上式成立，必须是第十章动力分析10.4解的稳定性10.4.2稳定性条件（10.57.21）（10