浅析各动力求解算法及其算法数值阻尼(人工阻尼)

唉呀妈呀是我呀

在谈数值阻尼之前，咱们先聊聊直接积分法的稳定性问题。

直接积分法是不对运动方程进行任何变换，直接运动方程进行积分求解。本质上讲，直接积分法是基于两个基本概念：

1.是将在求解域0≤t≤T内任何时刻t都应满足运动方程的要求，代之以仅在一定条件下近似地满足运动方程，例如仅在相隔△t的离散时间点0，△t，2△t，……上满足运动方程。

2.是在一定数目的△t区域内，假设位移D、速度V和加速度A的近似函数形式。

即在假定已知时刻t=0的位移D、速度V，并将时间求解域 0 ~ T等分成n个时间段0，△t，2△t，…，n△t（△t=T/n），且已求得时刻0，△t，2△t，…和t的解，现在要求t+△t时刻的解。建立了求解t+△t时刻解的算法后，即可类似地求得所有离散时刻的解，因此该方法也称逐步积分法。

逐步积分法是一种近似的求解方法，每一步积分计算都会带来误差。解的误差的来源主要有两大类：

第一类是在计算速度和加速度的近似表达式中略去了高阶小量，称为截断误差，这类误差是由方法本身决定的，一般可以作出估计，它随步长△t增大而增大。

第二类是由于计算机的字长总是有限位，超过其位数的数必须四舍五入，这类误差称为舍入误差。尽管在每一步中舍入误差并不大，但这类误差在求解过程中可能会被不断放大，以致使计算结果完全失真，这就是直接积分方法的稳定性问题。

如果无论△t取多大，给定任意初始条件，积分结果都不会无界增大，则称此方法是无条件稳定的。

反之，如果△t必须小于某个临界值时，积分结果才不会无界增大；则称此方法是条件稳定的。

在建立直接积分法格式时，我们假定已知时刻0，△t，2△t，…，t-△t和t的解，要求时刻t+△t的解，因此直接积分法利用时刻0，△t，2△t，…，t-△t和t的解递推求解t+△t时刻的解，即直接积分法的格式可以写成递推的形式

其中 Xt+Δt 和 Xt 为含有位移、速度等结果的向量，矩阵D称为递推矩阵或放大矩阵。

上式右端的第二项与载荷有关，在进行稳定性分析时可令L=0，此时给定初始条件X0后，时刻 t+△t = (n+1)△t 的解 Xt+Δt 可以写为

对矩阵D进行谱分解可得

假定矩阵未有重根(重根情况自行讨论)，则可得到

那么下式即为该递推矩阵的谱半径：

因此可以看出，若是谱半径≤1时，矩阵是有界的。综上所述，直接积分法的稳定性可归结为满足下式即为稳定：

如果满足以上稳定准则，当n→∞时，矩阵D^(n+1)有界。因此如果谱半径ρ(D)<1，当n→∞时，矩阵D^(n+1)→0，说明该积分方法存在数值阻尼（也称为人工阻尼）。ρ(D)越小，当n→∞时矩阵D^(n+1)趋于零的速度越快，表明该积分方法的数值阻尼越大。

因此谱半径除了度量算法的稳定性外，也度量了算法的数值阻尼。

对于一般结构动力学问题，系统的响应主要受低阶振型控制，高阶振型的贡献很小。

另外，由于受离散化的影响，有限元法或有限差分法对系统的高阶振型的近似程度很差。如果在直接积分中不能有效地滤除这些虚假的高阶分量，将会降低结果的精度。

接下来讨论下各算法的数值阻尼和用法的关键问题

01、【中心差分法】

【中心差分法】由于中心差分法所需要的时间步长比较短，实质上会让该算法的谱半径的模长等于1，也就是该算法并不能调整数值阻尼。

顺便提一句关于中心差分法的使用关键问题：对于n自由度系统，利用直接积分法对其运动方程进行积分，实质上等效于采用相同的时间步长 △t 同时对 n 个解耦后的微分方程进行积分。此时△t的选择应保证对每个微分方程的积分都是稳定的。

由于中心差分法的时间步长 △t 必须小于临界步长△tcr，因此它是条件稳定的。在用有限元法进行动力分析时，系统自由度数很高，结构系统的有效最小周期比较小，因此△t必须选得很小。另外在中心差分法中要求质量矩阵的全部对角元素都大于零，如果有对角元素等于零或接近于零，△tcr也将等于零或接近于零，中心差分法不能使用。

因此在用有限元法进行动力分析时，如使用中心差分法进行积分，则不宜使个别单元的尺寸比其他单元的尺寸小很多，否则会减小临界时间步长，因而大幅度增加计算量。

02、【Newmark-β法】

【Newmark-β法】Newmark-β法在适当的计算参数β、γ的选取下，是可以调整数值阻尼大小，甚至不要数值阻尼(如β=0，γ=1/2，则退化为中心差分法)。

Newmark-β法使用的关键问题：有图在导出Newmark积分格式时，使用的是t+△t时刻的运动方程，因此Newmark法是隐式方法。总所周知，Newmark法无条件稳定应满足：

如果不满足上述条件，要得到稳定的解，时间步长△t必须小于临界时间步长。

Newmark-β法的谱半径为：

当取为β=0，γ=1/2，此时Newmark法无数值阻尼，并具有二阶精度。只有当 γ>1/2 时，Newmark法才具有数值阻尼，其数值阻尼随着△t的增大而光滑地增大。但是，通过比较Newmark法的放大矩阵与系统的精确放大矩阵可知，此时Newmark法只有一阶精度(截断误差)，若是 Δt 较长的情况，即便计算也能趋于稳定，但结果误差也较大。为了使数值阻尼在高频段最大以有效地滤除系统虚假的高频响应，应使β取最小值，从而使ρ取最小值。

03、【Houbolt法】

【Houbolt法】Houbolt法是隐式计算方法，它的数值阻尼随 △t/T 的增大而很快增大，因此解的震荡衰减很快。当△t→∞时，Houbolt法的谱半径趋于0。

Houbolt法的使用关键问题，在Houbolt法中，不论△t取何值，在计算中Houbolt法的舍入误差都不会越来越大， Houbolt法是无条件稳定的。当M=C=0时该方法的递推方程可化为静力方程，因此Houbolt法可用来求解载荷随时间变化时的静态解。但是在中心差分法中当M=C=0时，有效质量阵M=0，因此中心差分法不能用来求解静态解，好在Houbolt法和中心差分法都是二阶精度。

04、【HHT-α法】

【HHT-α法】这个方法像Newmark和所有隐式方案一样，此方法的无条件稳定性适用于线性问题，HHT法的数值阻尼大小由α决定。

当 α = 0.0 时对应于Newmark方法，当满足以下式子时，HHT法对高频振动抑制最大。

其中α可以取0~0.33，且越大时，数值阻尼越大。
（注意：此处表达和OpenSees的α值表达不同，根据OpenSees文档中，α=1.0 对应Newmark法，和此处表达的内容是一样的。且对应OpenSees的HHT法为α建议取值再在0.67和1.0之间，且较小的α数值阻尼越大）

但是HHT方法中，即便是 Δt 较长的情况，计算也能趋于稳定，结果误差也不会很大。