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

        hello~

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

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

    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时刻解的算法后，即可类似地求得所有离散时刻的解，因此该方法也称逐步积分法。

    详情可以看下历史文章【JY】结构动力学初步-单质点结构的瞬态动力学分析。

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

    第一类是在计算速度和加速度的近似表达式中略去了高阶小量，称为截断误差，这类误差是由方法本身决定的，一般可以作出估计，它随步长△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)趋于零的速度越快，表明该积分方法的数值阻尼越大。

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

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

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

    举个例子，在【JY】结构动力学初步-单质点结构的瞬态动力学分析，这篇文章中，大型有限元软件Ansys和Opensees的计算结果出现了较多的虚假高阶分量，这篇文章做一个补充说明。

    可以看出Ansys和OpenSees等通用有限元软件通过计算时，对于绝对加速度计算时会产生虚假高阶高频分量（相对加速度影响不大，因为绝对加速度=相对加速度+地面加速度，会导致误差显著化。）

    因此一个好的直接积分方法应在高频段具有一定的可控的数值阻尼，以有效地滤除虚假的高频振型对系统响应的影响，同时在低频段的数值阻尼应尽可能小，以保证结果的精度。除了能滤除虚假的高阶振型的影响外，数值阻尼还有助于非线性问题迭代求解的收敛性，并且也有助于接触等具有约束的问题的求解。可见，在低频段谱半径 ρ(D) 应尽可能接近于1，而在高频段 ρ(D) 应逐步光滑地减小，并趋于某个给定值。

    Wood建议(Wood WL. Practical time-stepping schemes. Oxford: Clarendon,1990)，当△t/ T(结构周期)趋于无穷大时，ρ(D)趋于0.5~0.8较为合适。

（采用Hilber-Hughes-Taylor方法，添加数值阻尼）

    注意：结构在进行动力计算时，数值阻尼虽然能有效滤除虚假的高频振型，但同时会增加结构计算误差，使得结构输入输出总能量不同，多出数值阻尼产生的计算误差能量。因此，数值阻尼不宜过大，且要施加可控有效的数值阻尼（特别是单自由度或者等效单自由度结构，如：隔震结构），也可为抑制高频段的振动，增加刚度阻尼(详见：【JY】结构瑞利阻尼与经济订货模型)，该方法不会使得结构输入输出总能量不同，即结构能量是平衡的，但是计算手段视模型而定，不一定可行明显有效。

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

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

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

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

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

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

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

Newmark-β法的谱半径为：

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

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

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

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

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

    其中α可以取0~0.33，且越大时，数值阻尼越大。

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

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

建源工程侠，下期再见！

~路过别错过~