1、第四章：动力学问题的有限元方法,4.1动力学问题有限元方法的求解步骤,二维弹性动力学的基本方程是：,平衡方程：,几何方程：,本构方程：,边界条件：,初始条件：,结构在随时间变化的荷载作用下，位移、应变、应力等都是时间的函数。,建立静力学方程是从最小势能原理出发，而动力学方程是从Lagrange方程入手。,动力分析有限元方法的基本步骤：,1、连续区域的离散化：在动力分析中，因为引入了时间坐标，我们所处理的是（x,y,z,t)四维问题。在有限元分析中一般只对空间区域进行离散。,2、构造插值函数：因为只对空间离散，故单元内插值函数可表示为：,(4.1.1),或：,(4.1.2),其中：,注意：现在的
2、结点参数 是时间的函数。,3.形成单元、质量、刚度、阻尼、荷载矩阵,由位移对时间的微商是速度：,(4.1.3),根据lagrange方程建立动力学方程：,(4.1.4),式中,，,分别是系统的动能、势能、广义力,广义坐标和广义速度。,单元的动能和势能可分别表示如下：,（4.1.5),(4.1.6),（4.1.7),将（3.1.3）代入（3.1.5-7） 有,（单元质量矩阵）,（4.1.8),（单元刚度阵）,（单元荷载向量）,（单元阻尼矩阵）,（4.1.9),（4.1.10),这里 是整体结点位移向量，不同于前面的位移场 ，虽然记号相同。,4.形成整体求解方程,系统总动能为各单元动能之和：,系统
3、总势能为各单元势能之和：,系统总阻尼为各单元阻尼之和：,其中,其中,其中,这里整体刚度、质量、阻尼、荷载矩阵的形成方法和静力学刚度矩阵形成方法相同。,（总质量矩阵）,（总刚矩阵）,（总荷载）,（总阻尼阵）,将上述三式代入Lagrange方程有：,(4.1.11),这是一个二阶常微分方程组，若忽略阻尼的影响，方程可化简为：,(4.1.12),如果上式右端为“零”，则表示系统的自由振动。,（5）求解动力学方程（4.1.11).,（6）计算结构的应变和应力，一旦求得结点位移，便可由几何方程求得应变，进而求得应力等物理量。,4.2质量矩阵和阻尼矩阵,1、协调质量矩阵和集中质量矩阵,在上节中定义的单元质
4、量矩阵 称为协调质量矩阵或一致质量矩阵。这是因为导出它时，和导出刚度短阵所根据的原理(Gaterkin方法)及所采用位移插值函数是一致的同时质量分布也是按照实际分布情况考虑的。此外，在有限元法中还经常采用所谓集中(或团聚)质量矩阵。它假定单元的质量集中在结点上，这样得到的质量矩阵是对角线矩阵。,考虑平面应力应变单元，单元形式采用三结点三角形单元。,协调质量矩阵：,在三结点三角形单元中，位移的插值函数为：,（4.2.1),其中,是单元面积。,则有：,由三结点三角形单元的插值函数同其面积坐标相同的性质和面积坐标的积分公式：,容易得到三结点三角形单元的单元协调质量矩阵为：,其中,是单元的质量。,（4
5、.2.2),集中质量矩阵：,单元的每个结点上集中三分之一的质量，这样就得到单元的集中质量矩阵。,可以由单元的动能来推导单元的集中质量矩阵。,由单元的动能,写成矩阵的形式：,（集中质量矩阵）,即三结点三角形单元的集中质量为：,采用集中质量矩阵遇到的困难是对于高次单元如何将单元的质量分配到各个结点上，不像三结点三角形单元那样明显而简单，可能有多种选择，不易把握。另外，如果单元位移是协调的，同时单元刚度矩阵的积分也是精确的，则采用协调质量矩阵时，求得结构的频率将代表真实频率的上限，这点对设计工作是有意义的。,在实际分析中，协调质量短阵和集中质量矩阵都有应用，一般情况下，两者给出的结果也相差不多。质量
6、矩阵积分表达式的被积函数是插值函数的平方项，而刚度矩阵则是其导数的平方项，因此在相同精度要求条件下，质量短阵可用较低阶的插值函数，而集中质量矩阵从实质上看，正是这样一种替换方案。替换的好处是使计算得到简化，特别是采用直接积分的显式方案求解运动方程时，如果阻尼矩阵也采用对角短阵，可以省去等效刚度短阵的分解步骤，这点在非线性分析中特有更明显的意义。,2.振型阻尼矩阵,基于和协调质量矩阵的定义同样理由称,为协调阻尼矩阵。,在以后的讨论中，将知道系统的固有振型对于 和 是具有正交性的，因此固有振型对于比例于 和 的阻尼矩阵 也是具有正交性的。所以这种阻尼短阵称为比例阻尼或振型阻尼。今后还知道，利用系统
7、的振型矩阵对运动方程进行坐标变换时，振型阻尼短阵经变换后和质量短阵及刚度矩阵的情况相同，将是对角矩阵阵。这样一来，经变换后运动方程的各个自由度之间将是互不耦合的，因此每个方程可以独立地求解，这将对计算带来很大方便。,通常允许将结构的实际阻尼矩阵简化为 的线性组合即,其中 是不依赖于频率的常数，这种振型阻尼称为Rayleigh阻尼。,4.3直接积分法,直接积分是指在积分运动方程之前不进行方程形式的变换，而直接进行逐步数值积分。通常的直接积分法是基于两个概念，一是将在求解域 内的任何时刻 都应满足运动方程的要求，代之以仅在一定条件下近似地满足运动方程，例如可以仅在相隔 的离散的时间点满足运动方程。
8、二是在一定数目的 区域内，假设位移 ，速度 ，加速度 的函数形式。,在以下的讨论中，假定时间 的位移 ，速度 ，加速度 已知。并假定时间求解域 被等分为 个时间间隔 。 在讨论具体算法时，假定 时刻的解已经求得，计算的目的在于求 时刻的解。,1、中心差分法,将速度和加速度用位移的差分格式表示为：,（4.3.1),(4.3.2),在时刻 的位移解答 可由时刻 的运动方程应得到满足而建立，即：,(4.3.3),将（4.3.1)、(4.3.2)代入（4.3.3)有：,整理得：,(4.3.4),如果已经求得 、 ，则由上式可进一步解出 。,需要指出得是：该算法有一个起步问题。因为当 时，除了初始条件已
9、知的 而外，还需要知道 。,由,(4.3.5),削去 得,即,(4.3.6),而 可由 得到。,利用中心差分法求解运动方程得步骤为：,1.初始计算,（1）形成刚度矩阵 ，质量矩阵 和阻尼矩阵 。,（2）给定 ，并计算出 。,（3）选择时间步长 ， 并计算积分常数 、,。,（4）计算,（5）形成有效质量矩阵：,（6）三角分解,2.对于每一个时间步长：,（1）计算时间 的有效载荷,（2）求解时间 的位移,（3）如果需要，计算时间 的加速度和速度。,2.Newmark方法,Newmark积分方法实质上是线性加速度法的一种推广。它采用下列假设,(4.3.7),(4.3.8),在 时刻满足,(4.3.9
10、),由（4.3.8)得到：,(4.3.10),将其代入（4.3.7)有,（4.3.11),将（4.3.11、4.3.10)代入(4.3.9)有,Newmark方法具体算法：,1、初始计算,（1）形成刚度矩阵 ，质量矩阵 和阻尼矩阵 。,（2）给定 ，并计算出 。,（3）选择时间步长 ，参数 ，和 ，并计算积分常数：,（4）形成有效刚度阵,（6）三角分解,2 .对每一个时间步长：,（1）计算 时刻的荷载。,（2）求解时间 的位移,（3）求解时间 的加速度和速度,说明：Newmark方法是隐式方法，其稳定性与 的选取无关。,采用振型叠加法求解运动方程可分为以下三个主要步骤：,1、将运动方程转换到正
11、则振型坐标系,（1）求解系统的固有频率和固有振型,不考虑阻尼影响的系统自由振动方程：,4.4振型叠加法,它的解可以假设为以下形式,（4.4.1),（4.4.2),将（4.4.2)代入（4.4.1)得到广义特征值问题:,(4.4.3),其中,求解上式广义特征值问题,可以得到n个特征值和特征向量,分别记为:,称为系统的固有频率, 称为系统的固有振型.,固有频率和固有振型具有性质:,利用它们,原特征值问题可以表示为:,(4.4.4),(2)位移基向量的变换,因为 可以作为一个n 维空间的基底, 所以 可以按 展开:,其中,(4.4.5),将(4.4.5)代入(4.1.11)并左乘 有:,注意到特征向
12、量的正交性质有:,相应的初始条件转换成:,(4.4.6),两端左乘 有,(4.4.7),在(3.4.7)中如果如果阻尼是振型阻尼,则有 的正交性可得,写成矩阵的形式:,其中 是第 i 阶振型阻尼比.,在此种情况下, (3.4.6) 式为 n 个相互不耦合的二阶常微分方程,(4.4.8),上式中的每一个方程都相当于一个单自由度系统的振动方程, 可以比较方便的求解.式中 是载荷向量 在振型 上的投影. 若 是按一定的空间分布模式而随时间变化的,即:,则有:,上式中 表示空间坐标, 表示 在 上的投影, 为一个常数. 如 和 正交,则 ,从而得到 , .这表明结构响应中不包含 的成分. 亦即 不能激
13、起与 正交的振型 .需要求解的单自由度方程数也因之减少。,2、求解单自由度振动系统方程,单自由度系统的振动方程(3.4.8)的求解，在一般情况下可采用上节讨论的直接积分方法。但在振动分析中常常采用杜哈美(Duhamel)积分，又称为叠加积分。这个方法的基本思想是将任意激振力 分解为一系列微冲量的连续作用，分别求出系统对每个微冲量的响应，然后根据线性系统的叠加原理，将它们叠加起来，得到系统对任意激振的响应。杜哈美积分的结果是,(4.4.9),其中 是由初始条件决定的常数。,强迫振动项,自由振动项,当阻尼很小，即 时， ，这时杜哈美积分的结果是：,(4.4.10),3、振型叠加得到系统的响应,在得
14、到每个振型的响应以后，按(3.4.5)式将它们叠加起来，就得到系统的响应，亦即每个结点的位移值,振型叠加法的性质和特点：,对于n个单自由度系统运动方程的积分，比对联立方程组的直接积分节省计算费用。另外，通常只要对非耦合运动方程中的一小部分进行积分。例如只要得到对应于前 n 个特征解的响应，就能很好地近似系统的实际响应。,应该指出的一点是，如果在振型叠加法中，对于 n 个单自由度系统的运动方程都进行积分，且采用和直接积分法相同的积分方案和时间步长，则最后通过振型叠加得到的结果 和直接积分法得到的结果在积分方案的误差和计算机舍入误差的范围内将是一致的。,此外，对于非线性系统通常必须采用直接积分法。
15、因为此时 ，这样一来系统的特征解也将是随时间变化的，因此无法利用振型叠加法。,相应的特征值问题为,5.5缩减系统自由度的方法,在有限元动力分析中，发展提高计算效率、降低费用的数值方法是很有意义的。减缩系统自由度数日是广泛采用的方法之一。以下扼要地讨论两种减缩自由度的方法：主从自由度法。,1、主从自由度法,在主从自由度方法中将位移向量 划分为 和 两部分。并假定 按照一种确定的方法依赖于 。 称为主自由度， 称为从自由度。,用 表示 ：,(4.5.1),则有：,（4.5.2）,以无阻尼的自由振动方程为例：,（4.5.3),将(3.5.2)代入（3.5.3)有：,将上式左乘 有：,（4.5.4),其中：,下面来确定（3.5.1)中的 与 的关系矩阵 。,将 按静力方式施加于同一结构且不受其他荷载，由在结构内引起的变形模式确定 和 之间的关系。,（4.5.5),由静力平衡方程：,（4.5.6),由上式第二式得到：,即有：,这就得到了 和 之间的关系：,（4.5.7),将上式代入（3.5.5)有：,
