动力学方程数值解法：数值算法的性能分析

weixin

　　算法数值计算结果如何评价，针对不同的结构动力响应计算问题应该如何选择更合适的算法等是非常重要的问题。这就需要深入研究算法的数值计算性能，也就是算法的计算精度、稳定性等。

　　对线性结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼。

　　算法用于结构动力学方程的有限差分表示为：

　　以下算法的性能分析，均将算法用于这个方程。分别在相邻的不同时刻应用算法可得如下一般形式：

　　其中A为放大矩阵或称逼近算子，Lk为载荷逼近算子。

　　对自由振动情况有

　　显然计算的第n步的值与A直接有关。

　　例如，Newmak方法：

　　矩阵A的特征多项式为：

　　对Newmak方法有：

　　其中h为时间步长。

　　Wilson-θ方法，放大矩阵为：

　　放大矩阵A的特征多项式为：

　　其中A1，A2，A3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及矩阵的行列式，分别表达如下：

　　此外，在几个不同时刻应用数值算法，然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak方法，有：

　　1. 算法的稳定性分析
　　稳定性定义：设λi，其中i=1，2…m为放大矩阵A的特征值，则ρ=max∣λi∣定义为A的谱半径，若特征值互异，则ρ≤1的算法是稳定的，但若有重特征根，则要求ρ<1 。如果算法的稳定性要求对步长的选取有限制，称算法是有条件稳定的，反之为无条件稳定的。

　　判断方法：放大矩阵的谱半径小于等于1成立的充分条件是：

　　对3×3的放大矩阵：

　　上两式是关于算法自由参数ξ，Ω的不等式，由它可以判断算法是否无条件稳定，若不是，将给出稳定条件。

　　算法数值稳定性的物理解释：物理上，对一个无阻尼或者有阻尼自由振动系统，系统的能量随着时间不应该增加，有阻尼情况还应该减小。因此,一个数值方法的计算结果也不应该放大初始能量，如果经过若干步的数值计算以后，计算结果远比初始条件大，那就是数值算法本身计算是不稳定的。以下分析Newmak方法和Wilson-θ方法的稳定性问题。

　　Newmak方法：
　　将Newmak方法放大矩阵特征量代入稳定性分析表达式：

　　显然，当以下条件成立时，算法无条件稳定。

　　Wilson-θ方法：
　　对Wilson-θ方法有

　　上式很容易看出，其中第1、2、5不等式恒成立，对第3、4不等式若希望对任意的Ω均成立，则有：

　　求解上述不等式得

　　实际使用中通常选取θ=1.4。

　　算法的精度(相容性和收敛性)
　　直接积分算法的相容性、收敛性分析同样要使用其位移型的差分方程，用精确解代替近似解，得到局部截断误差表达式，用符号e(tk)表示

　　以最常用的线性三步法为例，局部截断误差用放大矩阵的特征量可表示为：

　　然后将下式在tk点进行泰勒展开，利用运动平衡方程化简即可。

　　算法精度定义：若局部截断误差表达式为步长的s阶小量，则称算法是s阶相容的。相容加稳定等于收敛，其相容的阶数就是算法的精度阶。收敛性的含义是当时间步长趋于零，算法的数值解趋于精确解。

　　下面来分析Newmak法的相容性和精度

　　其局部误差表达式得：

　　即将下式在tk时刻点泰勒展开。

　　并注意到在该时刻的运动方程有：

　　显然，当物理阻尼为零时，选择γ=1/2算法是二阶的。物理阻尼的存在，使算法精度降了一阶。

　　Newmak方法中有两个参数待定，每种特定的选取都是一个特定的算法，最常用的几个算法见表

　　常用的Newmak族直接积分算法

　　如果在一个时间步内需要求解一个隐式方程组，则称算法是隐式的，反之不需要求解方程，直接计算即可得到下一时刻的值，则称算法是显示的。

　　Newmak方法是隐式的，但对于中心差分方法，若质量矩阵和阻尼矩阵都是对角矩阵就可以显示地计算。显然显示方法计算量要小得多。Wilson-θ方法的精度，不难分析，无论是无阻尼还是有阻尼其精度都是2阶，隐式方法。

　　学习体会
　　通过这一学期计算结构动力的学习，使我对结构动力分析和数值主牌有了整体的了解。计算结构动力学论述了结构动力特性和响应的一般性原理和方法。对许多工程应用和计算机程序中常用的方法作了全面的介绍。

　　计算结构动力学是一门与各门其他力学课程相辅相成的学科，实际工程中的很多问题很复杂，对于计算出来的结果不一定能满足所有实际应用的要求。所以我们要重视这门学科发展争取在此基础上进一步提高。

　　本文摘录自百度文库《结构动力学方程常用数值解法》。