5．2．2 算法的稳定性分析 稳定性定义：设 i , i 1,2m 为放大矩阵 A 的特征 值， 则 maxi 定义为 A 的谱半径， 若特征值互异， 则 1 的算法是稳定的，但若有重特征根，则要求 1。 如果算法的稳定性要求对步长的选取有限制，称 算法是有条件稳定的，反之为无条件稳定的。 判断方法：放大矩阵的谱半径小于等于 1 成立的 充分条件是 1 2 A1 A2 0 1 2 A1 A2 0 1 A 0 2
其中 A1，A2，A3 为该矩阵的三个特征向量，分 别为矩阵的迹的一半、各阶主子式的和以及矩阵的 行列式，分别表达如下：
18 6 3 2 3 2 3 2 2 3 2 A1 2 2 2 ( 6) 4 2 3 2 3 6 2 2 18 12 A2 2 2 ( 6) 6 6 2 3 2 2 2 3 3 2 A3 2 2 ( 6)
C (a1 xt 2 xt a3 xt )
. ..
（2）求解 t t 时刻的位移
( LDL ) xt t Rt t
T

（3）计算在 t t 时刻的加速度、速度和位移
xt t a4 ( xt t xt ) a5 xt a6 xt
..
.
h 2 (6 3 2 2 2 2 ) 3 2 2 2 h(6 3 3 / 2) 6 2 3 3 2 2 6)
D ( 2 2 6)
放大矩阵 A 的特征多项式为：
det(A I ) 3 2 A12 A2 A3 0
例 5－1 分析 Newmak 方法、 Wilson- 方法的稳定性 解： 将 Newmak 方法放大矩阵特征量代入稳定性 分析表达式

1 2 ( ) 0 2
2
( ) (1 2 ) 1 0 2
2

1 6 2 0 , 2 显然，当 2 2 ， 6 (2 1) 0 2 2 算法无条件稳定。 12 ( 1 6 6 ) 0 对 Wilson- 方法有 (4 3 1 6 2 ) 2 24 12 0 4 (2 2 1 3 ) 0
1 1 A1 traceA ( A11 A22 ) A det A A A A A 11 22 12 21 , 2 2 2
对 Newmak 方法有：
v 1 2 [1 (2 1) ( ) ] 2 4 A1 D
1 2 [1 (2 2) ( ) ] 2 A2 D
..
xt t xt a7 ( xt t xt ) xt t xt txt a8 ( xt t xt )
5． 2 结构动力响应数值算法性能分析 算法数值计算结果如何评价，针对不同的结 构动力响应计算问题应该如何选择更合适的算法 等是非常重要的问题。这就需要深入研究算法的 数值计算性能，也就是算法的计算精度、稳定性 等。 对线性结构动力学问题，已经有证明对整个 多自由度的积分，等价于将模态分解后对单自由 度的积分的结果进行模态叠加，因此可以通过对 单自由度问题的分析，来说明算法的特性，其中 阻尼均假设为比例阻尼。
结 构 动 力 学
第五章 结构动力学数值算法
主讲教师：于开平
哈尔滨工业大学航天学院
5.1 结构动力学中常用数值方法
. .. M x C x Kx F (t ) x(0) x0 . x(0) v0
基本思想：首先给定待求时间长度 T， t [0, T ] ， 在其中取一系列离散点（i=0,1,…,n） ，我们不去求 x(t)，只求 x(ti ) 即可，即给出待求响应在各离散时 刻的近似值。两个离散时刻间隔称为步长，步长可 以相同（等步长） ，也可以不同（变步长） 。
T T
T T 2 x , hv 有的也 定义为 k k 或 xk , hvk , h ak
yk
yk 1 xk 1 , xk , xk m
对自由振动情况有
yn A y0
n
显然计算的第 n 步的值与 A 直接有关。 例如，Newmak 方法：
A A Ad
1 h At 2 h
计算数学领域，常微分数值算法常用的有两大类 1） 针对一阶微分方程数值积分法发展的欧拉法， 中点法，Rugge-kutta（龙格—库塔）方法。 2）直接基于二阶动力学方程发展的方法 对结构动力学问题的数值求解，常用的有两大类 1）模态迭加 2）直接积分 模态迭加方法， 比较常用， 但如下情况通常使用 直接积分方法（即求解之前不进行模态分析） i) 非比例阻尼，非线性情况。 ii) 有冲击作用，激起高频模态，力作用持续时间 较短，模态迭加计算量太大。
对 3 3 的放大矩阵
1 2 A1 2 A2 A3 0 3 2 A1 A2 3 A3 0 3 2 A1 A2 3 A3 0 1 2 A A A 0 1 2 3 1 A2 A3 (2 A1 A3 ) 0
通过变换将速度和加速度用位移表示，代入运动方 程，只剩 n+1 时刻位移一个未知数，得法 1 1 1） Newmark 平均加速度法， 2， 4 梯形公式 1 1 2） Newmark 线加速度法 2， 6 1 3） 中心差分法 2 ， 0
（2）求解 t t 时刻的位移
( LDL ) xt t Qt t
T
（3）计算 t t 时刻的加速度和速度
xt t a0 ( xt t xt ) a2 xt a3 xt
..
.
..
xt t xt a6 xt a7 xt t
.
.
..
上两式是关于算法自由参数 , 的不等式，由它可 以判断算法是否无条件稳定，若不是，将给出稳定 条件。
算法数值稳定性的物理解释： 物理上，对一个无阻尼或者有阻尼自由振动系 统，系统的能量随着时间不应该增加，有阻尼情况 还应该减小。 因此 , 一个数值方法的计算结果也不应该放大 初始能量，如果经过若干步的数值计算以后，计算 结果远比初始条件大，那就是数值算法本身计算是 不稳定的。
Newmark 类方法
1 2 1 3 x(tn 1 ) x(tn ) t x(tn ) t x(tn ) t x (tn ) O( t 4 ) 1） 2 6 . . .. 1 2 ... x(tn 1 ) x(tn ) t x(tn ) t x(tn ) O(t 3 ) 2） 2
.
1）可以直接略去高阶项 2）用变权来调节
1 xn 1 xn txn [( ) xn xn 1 ]t 2 2
xn1 xn [(1 ) xn xn1 ]t
然后假设在 tn 1 时刻近似满足运动方程
Mxn1 Cxn1 Kxn1 Fn1
2 x ( t ) x ( t ) t x ( t ) O ( t ) n n 3） n1 .. .. ...
x(tn1 ) x(tn ) x(tn ) O(t ) 由 3）得 代入 1） ，2）得 t
...
..
..
t 2 .. t 2 .. x(tn1 ) x(tn ) t x(tn ) x(tn ) x(tn 1 ) O(t 4 ) 3 6 . . t .. t .. x(tn 1 ) x(tn ) x(tn ) x(tn 1 ) O(t 3 ) 2 2
,
其中 h 为时间步长，
h, D 1 2
2
Wilson- 方法，放大矩阵为：
( 3 1) 2 6 1 A 3 D 2 6
h(6 2 3 2 ) ( 2 ( 2 3) 6) 6h
5.2.1 算法用于结构动力学方程的有限差分表示
2x 2 x f (t ) x
以下算法的性能分析， 均将算法用于这个方程。 分别在相邻的不同时刻应用算法可得如下一般形式
y k 1 Ayk Lk
A 为放大矩阵或称逼近算子， 为载荷逼近算子。
Lk
yk xk , xk 1 , xk m1
t
a4
a0
a5
a2

（4）形成等效刚度

K

K K a0 M a1C
（5）将等效刚度进行三角分解
K LDLT
2.对每一个时间步长 （1）计算 t t 时刻的等效载荷
. ..

Rt t Qt (Qt t Qt ) M (a0 xt a2 xt 2 xt )
纽马克法的解题步骤 初始值计算 （1）形成系统矩阵 K，M 和 C （2）定初始值 x0 ，
x0 ，
.
x0 。
..
（3）选择时间步长 t ，参数 、 （按上页选） 。 1 a0 a t 2 ， 1 t 并计算积分常数: t 1 a3 1 a4 1 a5 ( 2) 2 2 ， ，
此外，在几个不同时刻应用数值算法，然后将 方程中的速度和加速度项消去，可得数值算法关于 位移的差分方程，例如 Newmak 方法，有 1 2 2 (1 2 ) xn 1 2[1 (2 1) ( ) ] xn 2 4
1 2 [1 (2 2) ( ) ]xn 1 0 2
..
5.1.2 威尔逊- 法的解题步骤 1. 初始值计算 （1）形成系统矩阵 K，M 和 C （2）定初始值 x0 ， x 0 ， x 0 。 （3）选择时间步长，并计算积分常数
.
..