L o a d i n g . . . . . .
0 Applause for me

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

        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】结构瑞利阻尼与经济订货模型),该方法不会使得结构输入输出总能量不同,即结构能量是平衡的,但是计算手段视模型而定,不一定可行明显有效。

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

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

建源工程侠,下期再见!

往期精彩*点击即达

#性能分析

【JY】基于性能的抗震设计浅析(一)

【JY】基于性能的抗震设计浅析(二)

【JY】浅析消能附加阻尼比

【JY】近断层结构设计策略分析与讨论

#概念机理

【JY】基于Ramberg-Osgood本构模型的双线性计算分析

【JY】结构瑞利阻尼与经济订货模型

【JY】结构动力学初步-单质点结构的瞬态动力学分析

【JY】从一根悬臂梁说起

【JY】反应谱的详解与介绍

【JY】主成分分析与振型分解

【JY】浅谈结构多点激励之概念机理(上)

【JY】浅谈结构多点激励之分析方法(下)

【JY】板壳单元的分析详解

#其他

【JY】位移角还是有害位移角?

【JY】如何利用python来编写GUI?

~路过别错过~

欢迎关注微信公众号: 建源学堂