第三章 离散化结构动力方程的解法.doc

第三章 离散化结构动力方程的解法3.1 绪 言对于一个实际结构，由有限元法离散化处理后，应用瞬时最小势能原理可导出动力方程 （3.1）这里，、及分别表示加速度、速度、位移及所作用的外力矢量，他们都是与时间有关的。从数学的角度来看，式（3.1）是一个常系数的二阶线性常微分方程组，对于它的求解原则上并无困难。但是，由于、和的阶数非常高，使得式（3.1）的求解必须花费很大的代价，便促使人们去寻求一些效率高的近似计算方法。目前，用于求解式（3.1）的方法，大致可分为两大类。一是坐标变换法，它是对结构动力方程式（3.1），在求解之前，进行模态坐标变换，实际上就是一种Ritz变换，即把原物理空间的动力方程变换到模态空间中去求解。现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。通过这一变换，实现降阶，求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。显然，这种方法的计算精度，取决于所假设的模态。用Ritz矢量法求解的近似模态作为假设模态，可得到满足要求的精度。二是直接积分法，它是对式（3.1）在求解之前，不进行坐标变换，直接进行数值积分计算。这种方法的特点是对时域进行离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成，于是，式（3.1）就化为一个由位移组成的该离散时刻上的响应值，通常又称为逐步积分法。线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。主要有中央差分法，Houbolt方法，Wilson-法和Newmark方法等。3.2 模态（振型）迭加法设有n个自由度的系统，在外力的作用下，常常被激起较低阶的一部分模态（即振型），而绝大部分高阶模态被激起的分量很小，一般可忽略不计。例如，在地震载荷作用下，通常，只有最低的二阶，三阶模态起主要作用。所以，对于这样的一些问题，采用模态迭加法是有效的。设有式（3.1）的n阶动力方程，起主要作用的是其前q阶模态，通常取。按Ritz变换，则可将式（3.1）中的用前q个模态的线性组合来表示，即 其中，为结构的已知的保留主模态矩阵，而是维的模态基坐标矢量，它形成了一个q维的模态空间。它表示在中，各阶主模态所占有的成分的多少。 假定已用第二章所述的某一方法解出，再将式（3.2）代入（3.1），并左乘以，可得 式中显然，式（3.3）是一个q阶的微分方程组。由于，所以，它比式（3.1）的n阶就小的多了，实现了降阶，因而也就容易求解多了。 若展开上述的的表达式，根据主模态（主振型）关于的表达式，根据主模态的（主振型）关于M的正交性质，可知所以，是一个对角阵。同理可知也是一个对角阵。然而，在一般的情况下，是一个非对角阵，即在模态空间中，系统的的阻尼一般是耦合的。因此，式（3.3）是一个完全解耦的动力学方程。但是，它是一个已降阶的q阶的动力方程，可使用后面即将介绍的直接积分法求解。 当系统的阻尼为比例阻尼时，即可以表示为 则为对角阵。此外，若系统的阻尼是一般的的线性阻尼，并非比例阻尼，但是只要结构的固有频率不相等，而且不十分接近，则可用舍去阵中的非对角元来实现的对角阵，也不会引起太大的误差。在上述两种情况下，可以获得对于模态坐标的完全解耦的动力学方程。即式（3.3）是q个独立的方程，每个方程只包含一个未知量，相互之间不耦合。因而式（3.3）可按单自由度的动力学方程写为 或 其中。式（3.6）可用直接积分法计算，或用Duhamel积分求得其解为 式中，而，由初始条件 得出的 与决定。 由于有阻尼的存在，由初始条件所激发的振动，随时间的增长而衰减以致消失。因此，常可不计式（3.7）中的第二项，即是由初始条件激发的自由衰减振动。计算出后，便可利用式（3.2），计算出物理坐标的响应。整个计算步骤可归纳如下:第一步：根据结构的离散化模型，建立系统的以及，并进行结构的固有特性分析，即求解特征值问题 求出前阶特征对，（）第二步：形成模态阵，并建立模态基坐标下的动力方程其中，而。根据实验结果或经验数据确定各阶主振动中的比例阻尼。第三步：求解主模态基坐标的动力方程，有，其中，。第四步：进行坐标变换后，求得动力响应3.3模态假设法上节所述的模态迭加法，是用系统的真实主模态组成的模态矩阵，再对系统的物理坐标进行模态坐标变换，从而在主模态空间中得到降阶并解耦的动力学方程，这样来实现简化计算。而这里提出的假设模态法，则是用一组假设模态矩阵，对系统的物理坐标进行模态坐标转换，从而在模态空间中得到一组只降阶的动力学方程。若令假设模态矩阵为，而，进行坐标变换，即 (3.10)把它代入式（3.1），并左乘，则可得到降阶的动力学方程为 (3.11)其中， , 。它们分别对应于假设模态坐标的质量矩阵、阻尼矩阵、刚度矩阵与广义力列阵。因为矩阵中的各列都是假设模态，它们一般不具有正交性，所以，、和都不是对角阵。于是，方程(3.11)是不能解耦的方程组，但它却是比式(3.1)的阶数要低得多了。显然，对式(3.11)采用直接积分法求解，将比对式(3.1)求解要简便得多。这是假设模态法的优点。假设模态法的计算精度，很显然地是取决于假设模态阵中模态假设的好坏与质量。因此，应用假设模态法能否成功的关键在于确定出一个适宜的假设模态矩阵。在第五章中，我们介绍了几种构造假设模态的方法。实际上，Rayleigh-Ritz法可认为是一种假设模态法。它的作用在于降低方程的阶数，简化计算。它的基本思想是，事先假定出若干近似的特征矢量，然后按照这些特征矢量的最佳线性组合，而算得前若干阶特征值的近似值。显然，运用这种方法时，其计算精度与事先假定的特征矢量的近似程度和数量有关。按照Ritz变换的思想，找到了近似的特征矢量后，即有 (3.12)求解如下的广义特征值问题，即 (3.13)其中和为原结构离散化之刚度阵和质量阵，它们都是阶方阵。求解式(3.13)，得到个特征矢量，有再按照Ritz的变换，即式(3.12)，由特征矢量，可计算出矢量，即是 (3.14)现在用来表示此变换阵，它就是我们要构造的假设模态矩阵。3.4 中心差分法（显示法）现在开始讨论直接积分法，或称逐步积分法。前面讨论的模态迭加法，并非总是有效的。当刚度矩阵，或质量矩阵，或阻尼矩阵出现随时间变化时，或当外荷载激起的振型太多，需要计算的特征对太大时，就不宜于采用模态迭加法，在这些情况下，采用逐步积分法是适宜的。中心差分法就是其中的一种。这种方法的特点，是将动力方程在时间域上离散，化成对时间的差分格式，然后根据初始条件，利用直接积分法逐步求解出一系列时刻上的响应值。假定时，位移、速度和加速度分别为已知的，和。再将求解的时间区间划分为个等分，即。我们要建立的积分格式就是从已知的，的解来计算下一个时间步的解。在中心差分法中，按中心差分将速度和加速度矢量离散化为（3.15）（3.16）上面二式将时刻的速度和加速度用相邻时刻的位移表示了。考虑在时刻的动力方程，有（3.17）将式（3.15）和（3.16）代入式（3.17）中，得到（3.18）这样，上式就化为用相邻时刻的位移表示的代数方程组。由它可解出。又由于它是利用时刻的方程解得的，所以，它称为显示积分。并且，还注意到，在求解时，需要用，的值。于是，在计算开始时，即时，要计算的值、就需要的值，他是未知的，因此，必须有一个启动的处理，因而这种算法不是自起步的。由于和是已知的，所以，由时的式（3.15）和（3.16），可解得（3.19）使用中心差分法的逐步求解过程如下：A初始计算（1）形成刚度矩阵，质量矩阵和阻尼矩阵。（2）给定初始值和。（3）选择时间步长，并计算积分常数：，。（4）计算。（5）形成有效质量矩阵。（6）三角分解：。B对每个时间步计算（1）计算时刻的有效载荷。（2）求解时刻的位移。（3）如果需要计算时刻的速度和加速度应当指出，这种中央差分算法，左端的系数矩阵只与质量阵和阻尼阵有关，而与刚度阵无关。如果质量阵和阻尼阵是对角阵，那么在解方程时，就不需要对系数阵进行三角分解，即不需要解线性代数方程组，从第一步开始逐次直接求得各个时刻的值，这是中央差分格式就是一种显示的格式。此外，由于不求解代数方程组，也就不需要进行组集，它的右端项的形成也只须在单元一级水平上，由每个单元对有效载荷矢量的贡献迭加而成。因此，ADINA程序规定，在用中心差分法时，必须使用对角的质量阵和阻尼阵。从计算稳定性角度来看，中心差分法的缺点，在于它是条件稳定的，即当时间步长太大 时，积分是不稳定的。所以，对步长的限制是这里，是临界步长值，是有限元系统的最小周期。这样，当很小时，就限制了必须很小，所以求解所花的代价就很大。3.5 线性加速度法和Wilson-法线性加速度法和Wilson -法，都是属于逐步积分法。线性加速度法是假定在时间间隔内，即在步长时间内，加速度呈线性变化，其表达式为（3.20）其中，。但是，这个方法不是无条件稳定的，所以在应用上受到限制。70年代初期，Wilson推广了线性加速度法，他假定在此步长更大的时间区间内，加速度仍保持线性变化，经过证明，当时，这一方法是无条件稳定的，这就是Wilson -方法。这个方法的加速度表达式为(3.21)式中显然，对比式（3.20）和式（3.21）得知，线性加速度法是Wilson -法中，当时的一个特例。所以，我们只讨论Wilson -法就够了。在区间内，对式（3.21）进行积分，得到（3.22）和（3.23）令，由上二式，有（3.24）和 （3.25）从这二式，可将时刻的加速度和速度用位移来表示即（3.26）和(3.27)于是，在时刻的动力方程为(3.28)式中，将(3.26)和式(3.27)代入式(3.28)，就得到关于的方程为(3.29)记于是，式(3.29)可写为(3.30)求解方程(3.30)，则得到将求解得到的，代入(3.26)中，就得到。如在(3.21)中，取，并将式(3.26)代入，有(3.31)将(3.21)代入式(3.22)和(3.23)，并取，有(3.32)(3.33)用法逐步求解的过程如下：A 初始计算（1） 形成刚度矩阵，质量矩阵和阻尼矩阵。（2） 给出初始值，和。（3） 选择时间步长，取，并计算积分常数，, ， ， ， ， , , （4）形成有效刚度矩阵：（5）对作三角分解：B对每个时间步计算（1）计算时刻的有效载荷 （2）计算时刻的位移 （3）计算时刻的位移，速度和加速度 与中心差分法相比较，Wilson-法是隐式积分，即每计算一步，必须解一个线性代数方程组。当时，它是无条件稳定的。此外，这种算法是自起步的，时刻的位移，速度和加速度都可由时刻的变量表示，不需要特别的起动处理。 3.6 Newmark方法Newmark在1959年提出的逐步积分格式，故称为Newmark方法。它的基本假定是 (3.34) (3.35)其中和是按积分的精度和稳定性要求可以调整的参数。当, 时，它就是线性加速度法，所以，Newmark方法也可以理解为线性加速度法的一个小延伸。Newmark法最初提出作为无条件稳定的一种积分格式是常平均加速度法，即假定从到时刻，加速度不变，取为常数。此时，取, 。 常平均加速度法是应用得最广泛的逐步积分方法之一。研究表明，当，时，Newmark方法是无条件稳定的。从式和可得到， 用及、和表示的表达式，即有 （3.36）和 （3.37）考虑时刻的动力方程，有 （3.38）将式（3.36）和（3.37）代入（3.38），就得到关于的方程为 （3.39）其中求解方程（3.39），就可得到，然后，根据式（3.36）和式（3.37）可解出和。Newmark方法逐步求解的过程如下：A. 初步计算 （1）形成刚度矩阵，质量矩阵和阻尼矩阵。 （2）给定初始值，和 （3）选择时间步长，参数和，并计算积分常数： （4）形成有效刚度矩阵 （5）对作三角分解： B对每个时间步计算 （1）计算时刻的有效载荷 （2）求解时刻的加速度和速度我们注意到Wilson-法与Newmark法的计算关系式，在形式上是相同的，只是其中的系数取不同的值而已。因此，它们可用同一计算机程序来实现。3.7 Houbolt方法这个差分格式是利用、四个时刻上位移的三次插值多项式建立起来的。即假定 （3.40）和 （3.41）这里认为，和是已知的，而是未知的。考虑时刻的动力方程，有 （3.42）将式（3.40）和（3.41）代入式（3.42）中，就得到求解时刻的方程为 （3.43）由上式解得后，代入式（3.40）和式（3.41）中，便求得了和。这样，逐步求下去，便可求得任意时刻的动力响应值。应该注意到，这个差分格式不是自起步的，除了利用初始条件，和之外，尚需利用前述的任一种自起步的方法，求得、和后，再利用式（3.40）和式（3.41），即由上二式，可求 、。这样，利用，和，就可由开始