唉呀妈呀是我呀 发表于 2022-9-1 13:26

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

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

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

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 较长的情况,计算也能趋于稳定,结果误差也不会很大。

页: [1]
查看完整版本: 浅析各动力求解算法及其算法数值阻尼(人工阻尼)