动力学常微分方程的时间积分方法(精)/计算力学前沿丛书

第1章 α类时间积分方法 Newmark方法[1]是在工程中*先得到广泛应用的时间积分法，它包含了很多著名的格式，如中心差分法，梯形法则和 Fox-Goodwin方法[2]。这些格式均为具有二阶精度的辛几何格式，它们完全没有数值耗散，不会在系统中引入额外的数值阻尼。但值得说明的是，对于大型有限元系统，空间离散造成的虚假高频模态的参与可能会破坏数值解的准确性，甚至导致矩阵病态，无法得到收敛解。因此，数值阻尼在某些情况下对于给出稳态解或完成暂态全过程分析是十分必要的。*早提出的多步法，如 Houbolt方法[3]和 Park方法[4]，虽然拥有强烈的高频耗散能力，但它们的精度较差，且需要额外的启动程序，在实际问题中很少使用。 Wilson-θ方法[5]是*早的具有可控数值阻尼的时间积分法，它可以实现无条件稳定和二阶精度，但是 Goudreau 和 Taylor[6]指出 Wilson-θ方法关于初始位移有二次超调，关于初始速度有一次超调，这使得它在用于计算非零初始值问题时，在开始的几步可能会出现初始数据及其误差被放大的现象。此外， Newmark方法虽然具有耗散格式，但其精度为一阶并且数值阻尼不可控制。为了改善这一性能，HHT-α方法[7]，WBZ-α方法[8]，HP-θ方法[9]，广义α方法[10，11]和各类耗散算法[12]陆续被提出，这些方法均具有二阶精度，无条件稳定性，以及可控的高频耗散能力。其中，由于广义α方法具有更多的参数组合形式，在相同的高频耗散程度下可以实现更高的低频精度，在实际应用中比较受欢迎。 广义α方法沿用了 Newmark方法的假设，但它往平衡方程中引入了额外的参数，使得平衡方程在时间节点处不严格满足，而是在时间区间或步长内加权满足。本书作者通过引入额外变量的方式，提出了一种新的三参数方法[13]，它具有和广义α方法相同的谱特性，但由于它严格满足平衡方程，避免了载荷插值引入的误差，具有一定的加速度精度优势。采用类似的思想，本书作者还提出了一种四参数方法[14]，通过选取合适的参数值，它*高可以达到四阶精度，而且不具备数值耗散，适用于长期动力学仿真。 1.1 Newmark方法 结构动力学问题可以用如下常微分方程的初值问题来描述 (1.1.1) 其中，M 为质量矩阵，N 为结构内力，包含了阻尼力和弹性力等恢复力，R 为外部载荷，和分别表示位移，速度和加速度向量，t 表示时间，x0和 v0为已知的初始位移和速度。对于线性系统，该方程可进一步简化为 (1.1.2) 其中，C 和 K 分别为阻尼和刚度矩阵。一般来说，对结构线性动力学问题的求解，常用的有模态叠加法和时间积分法两大类。模态叠加法是建立在坐标变换基础上的解析方法，实际上是常微分方程由特解求通解的方法。对于复杂多自由度系统，求出全部模态或频率基本上是不可能的，对于实际问题也是没有必要的，因此只有对自由度较少的系统，才可能用全部的模态来叠加得到解析解。相比较来说，时间积分法是更加通用的求解技术，它不需要坐标变换，适用于任意激励或非线性情况，在工程实际中得到了广泛的应用。时间积分法的基本思想是：首先给定待求时间长度 T，即 t ∈[0， T]，在其中布置一系列离散时间点 ti (i =0，1， ， N)，使用差分格式由已知量给出待求变量在下一离散时刻的近似值，动力学平衡方程在时间点上得到满足或在相邻点用加权形式来满足。不同的差分格式对应不同的时间积分法，也给这些方法赋予了不一样的数值性能。 1.1.1算法格式 Newmark方法是结构动力学中应用*广的时间积分法，它采用的差分格式如下 (1.1.3) 其中，β和γ为两个自由参数，h 为时间步长，下标“k”表示状态变量在 tk 时刻的近似值。Newmark方法还用到了在时间节点 tk+1处的平衡方程 (1.1.4) 为启动计算程序，在给定初始位移和速度的情况下，初始加速度需由平衡方程得 (1.1.5) 每一步需要求解的方程为 (1.1.6) 由式(1.1.6)得到第(k+1)步的加速度之后，代入方程(1.1.3)，则可得到第(k+1)步的位移和速度。若系统非线性，则每一步需要求解一个非线性代数方程，这需要借助 Newton-Raphson 等迭代方法来实现。表1.1中给出了 Newmark方法用于线性系统的计算流程，对于大型有限元系统，它的计算量主要花费在矩阵分解运算上。当β=0时，从表1.1中可以看出，此时 Newmark方法不需要分解刚度矩阵，这对于质量矩阵和阻尼矩阵均为对角阵的情况十分有利，因此称β=0时为显式格式，其余情况则为隐式格式。 表1.1 Newmark方法用于线性系统的计算流程 A.初始准备 1.空间离散得到质量矩阵 M、阻尼矩阵 C 和刚度矩阵 K； 2.初始化 x0和 x˙0，计算初始加速度； 3.选取参数γ和β，以及时间步长 h； 4.建立常量矩阵，并将其三角分解。 B.第(k +1)步 1.计算有效载荷向量； 2.计算加速度； 3.计算位移； 4.计算速度。 1.1.2 数值性能 如何评价时间积分法的计算结果，针对不同的问题如何选择合适的算法，这些都需要评估时间积分法的数值性能，包括稳定性、精度、超调、数值耗散和弥散等。本章以 Newmark方法为例，给出精度和稳定性分析的一般流程。 目前来看，时间积分法的性能分析仍局限于线性动力学范畴。对于线性多自由度系统的积分，等价于将其模态分解后对单自由度系统积分的结果进行模态叠加。因此可以通过对单自由度问题的分析来说明算法的特性，即用于性能分析的模型方程为 (1.1.7) 其中，ξ、ω和 r 分别为阻尼率、固有频率和模态分解后的外激励。应用于该单自由度系统，Newmark方法的递推格式可以简化为 (1.1.8) 其中，为载荷作用向量，A 为 Jacobi 矩阵，它对算法性能起着决定性作用，其形式为. (1.1.9) 其中，τ=ωh。围绕着 Jacobi 矩阵 A，Bathe 和 Wilson[15]给出了精度阶数的定义和稳定性的判别方法，接下来以 Newmark方法为例一一进行介绍。 (1)稳定性 对于一个稳定的自由振动系统，系统能量随着时间不应该增加，有正阻尼情况还应该减小。因此，差分方法的计算结果也不应该放大初始能量。如果经过若干步的数值计算后，计算结果远比初始条件大，甚至随着时间的增长发散到无穷，那就说明数值算法本身是不稳定的。令 r =0，由方程(1.1.8)可以得到 (1.1.10) 当k→∞时，Xk是有界的当且仅当 (1.1.11) 式(1.1.11)是算法稳定的充要条件。若稳定性要求对时间步长的选取有限制，则称该时间积分法是条件稳定的，反之为无条件稳定。 针对Newmark方法，矩阵A的本征多项式可写为 (1.1.12) 其中，A1、A2和 A3分别为矩阵A的迹、二阶顺序主子式之和以及矩阵的秩。由方程(1.1.9)可以得到 (1.1.13) (1.1.14) (1.1.15) 从方程(1.1.15)可知，矩阵 A 有一个零本征根。根据 Routh-Hurwitz稳定性判据[16，17]，其余两个本征根的模小于等于1的充要条件为 (1.1.16) 从而可以得到Newmark方法的无条件稳定条件为 (1.1.17) 而条件稳定性要求 (1.1.18) 其中，τcr 称为稳定极限。对于条件稳定算法，系统的*高频率决定了步长的可取范围，因此，大型有限元系统的可取步长可能很小，这使得无条件稳定算法在结构动力学问题分析中通常更受欢迎。稳定性要求是一个算法能够得到使用的前提条件，不稳定的算法无法给出合理的数值解。 (2)精度 根据 Lax 定理[18]，收敛的算法除了要满足稳定性要求外，还应满足相容性条件，即算法应至少具有一阶精度。令 r =0，在Newmark 算法格式(1.1.8)中利用平衡方程消去速度和加速度项可以得到 (1.1.19) 它的局部截断误差定义为 (1.1.20) 若σ= O(hl)，则称该方法具有 l 阶精度；若1.1，则称该差分格式(1.1.8)与微分方程(1.1.7)相容。将式(1.1.20)在 tk 处级数展开，并将方程(1.1.13)～式(1.1.15)代入，可以得到 (1.1.21) 因此，Newmark方法至少具有一阶精度，若γ=1/2，它可以达到二阶精度。特别地，对于无阻尼系统，即ξ=0，若γ=1/2且β=1/12，Newmark方法可以实现四阶精度，这就是著名的 Fox-Goodwin方法[2]。 当时间步长趋于零时，收敛算法给出的数值解会逐渐趋近于解析解。精度阶数越高，数值解与解析解之间的差别越小。但是超过二阶精度的多步和单步算法往往无法实现无条件稳定[19]，因此常用的时间积分法大都为二阶精度。 Newmark方法包含了几种著名的格式，列举如下： 1.中心差分法(γ=1/2，β=0)(Central Difference Method) 中心差分法是一种条件稳定的显式辛几何方法，稳定极限为τcr =2。当质量和阻尼矩阵为对角阵时，中心差分法仅需做向量运算，无须进行矩阵分解，效率较高，广泛用于求解波传播、冲击和非线性等问题。 2.梯形法则(平均加速度方法，γ=1/2，β=1/4)(Trapezoidal Rule) 梯形法则是无条件稳定的隐式方法，对时间步长的选取没有限制。当它用于线性无阻尼系统时，梯形法则可以严格保守系统的能量，与欧拉中点辛差分格式等价，是一种隐式辛算法。 3. Fox-Goodwin方法(γ=1/2，β=1/12) Fox-Goodwin方法是一种条件稳定的隐式方法，稳定极限为。Fox- Goodwin方法的优势在于它用于线性无阻尼系统时精度较高，可以达到四阶精度。 值得注意的是，当γ=1/2时，对于无阻尼系统，Newmark方法的谱半径在稳定区间内保持为1，也就是说，Newmark方法的二阶格式完全不具备数值阻尼，这使得它在求解刚性问题和一些非线性系统时遇到了困难。 1.2 广义α方法 为了提高 Newmark方法的阻尼性能，人们陆续发展了一些具有可控数值耗散性能的二阶方法，包括 Wilson-θ方法[5]、HHT-α方法[7]、WBZ-α方法[8]、HP-θ方法[9]、广义α方法[10，11]等，本章以广义α方法为代表进行介绍。 1.2.1 算法格式 广义α方法沿用了 Newmark方法的差分格式