?n?1?x?n?[(1??)????xxn??xn?1]?t然后假设在tn?1时刻近似满足运动方程
通过变换将速度和加速度用位移表示,代入运动方程,只剩n+1时刻位移一个未知数,得
参数不同选取包含着三个经典算法 (1)Newmark平均加速度法
??M?xn?1?Cxn?1?Kxn?1?Fn?1Kxn?1?Qn?1???11??2,4
(2)Newmark线加速度法
??11??2,6
(3)中心差分法
??12,??0
Newmark法的一般步骤: 1初始值计算
(1)形成系统矩阵K,M和C (2)定初始值x0,x0和x0
(3)选择时间步长?t,参数?、?。并计算积分常数:
a0?1?1??t?a?a??1a??1a?(?2)45132a??t(1??),a7???t ??t,?,2???t 2?,,6?...(4)形成等效刚度矩阵K
K?K?a0M?a1C(5)K???
矩阵进行三角分解
TK?LDL 2对每一时间步
(1)计算t??t时刻的等效载荷
Qt??t?Qt?M(a6xt?a2xt?a3xt)?C(a1xt?a4xt?a5xt)(2)求解t??t时刻的位移
?......
(LDL)xt??t?Qt??t(3)计算t??t时刻的加速度和速度
T?
xt??t?a0(xt??t?xt)?a2xt?a3xt.....xt??t?xt?a6xt?a7xt??t
Wilson-?法的一般步骤: 1初始值计算
(1)形成系统矩阵K,M和C (2)定初始值
......x0,
x0.,
x0。
..(3)选择时间步长,并计算积分常数
6a?a3??t3?ta0?a1?a3?a4?0a5?2a6?1?a7?2a2?2a1,(??t),??1.4,??t,2,?,?,?,2,
?t2a8?6。
(4)形成等效刚度K
?K?K?a0M?a1C(5)将等效刚度进行三角分解
TK?LDL ??
2对每一个时间步长
(1)计算t??t时刻的等效载荷
Rt???t?Qt??(Qt??t?Qt)?M(a0xt?a2xt?2xt)?C(a1xt?2xt?a3xt)
(2)求解t???t时刻的位移
?......
(LDL)xt???t?Rt???t..T?
...(3)计算在t??t时刻的加速度、速度和位移
xt??t?a4(xt???t?xt)?a5xt?a6xt?t??t?x?t?a7(??t??t???t)xxx?t?a8(??t??t???t)xt??t?xt??txxx
三 结构动力响应数值算法性能分析
算法数值计算结果如何评价,针对不同的结构动力响应计算问题应该如何选择更合适的算法等是非常重要的问题。这就需要深入研究算法的数值计算性能,也就是算法的计算精度、稳定性等。
对线性结构动力学问题,已经有证明对整个多自由度的积分,等价于将模态分解后对单自由度的积分的结果进行模态叠加,因此可以通过对单自由度问题的分析,来说明算法的特性,其中阻尼均假设为比例阻尼。 算法用于结构动力学方程的有限差分表示为:
???2??x???2x?f(t) x以下算法的性能分析,均将算法用于这个方程。分别在相邻的不同时刻应用算法可得如下一般形式:
yk?1?Ayk?Lk
其中A为放大矩阵或称逼近算子,Lk为载荷逼近算子。
yk??xk,xk?1,?xk?m?1?对自由振动情况有
TTyk?1??xk?1,xk,?xk?m?yn?Any0
显然计算的第n步的值与A直接有关。 例如,Newmak方法:
A?At?1Ad
?h22h[1?h(1?2?)??]?1?(1?2?)?222??1?h??2h????Ad?2At???1?2h(1??)????2?21?2h?????h(1??)???h??? ,
矩阵A的特征多项式为
?det(A??I)??2?2A1??A2?0
A1?11traceA?(A11?A22)22,
A2?detA?A11A22?A12A21对Newmak方法有:
[1?(2??1)???(??A1?Dv1?)?2]24,
[1?(2??2)???(????A2?其中 h为时间步长,
1)?2]2
D???h,D?1?2??????2
Wilson-?方法,放大矩阵为:
?(??1)??6?1?A???3??D??6?2?h2(6???3?2?2??2?2)32?h(6?????)?(?(??3)?6)6h???22
232??3222h(6?????3?3??/2)??6???2?3?3?2?2?6)?
D?(?2?2?6)?
放大矩阵A的特征多项式为:
?det(A??I)??3?2A1?2?A2??A3?0
其中A1,A2,A3为该矩阵的三个特征向量,分别为矩阵的迹的一半、各阶主子式的和以
及矩阵的行列式,分别表达如下:
此外,在几个不同时刻应用数值算法,然后将方程中的速度和加速度项消去,可得数值算法关于位移的差分方程,例如Newmak方法,有
?18??6?3?2?3??2?3?2?2?3?2??A1?222?(???6)??4?2?3?2?3?6?2?2?18??12??A2?22?(???6)???6?6???2?3?2?2??2?3?3?2??A3?22?(???6)??(1?2??????2)xn?1?2[1?(2??1)???(???2?1)?2]xn4
?[1?(2??2)???(????1)?2]xn?1?02
3.1算法的稳定性分析
稳定性定义:设
?i,i?1,2?m为放大矩阵A的特征值,则??max?i定义为A的谱半径,
若特征值互异,则??1的算法是稳定的,但若有重特征根,则要求??1。如果算法的稳定性要求对步长的选取有限制,称算法是有条件稳定的,反之为无条件稳定的。
判断方法:放大矩阵的谱半径小于等于1成立的充分条件是
1?2A?1?A2?0?1?2A?1?A2?0?1?A2?0?
对3?3的放大矩阵
1?2A?1?2A2?A3?0?1?A2?3A3?0?3?2A??3?2A1?A2?3A3?0?1?2A1?A2?A3?0??1?A2?A3(2A1?A3)?0?
上两式是关于算法自由参数?,?的不等式,由它可以判断算法是否无条件稳定,若不是,将给出稳定条件。
算法数值稳定性的物理解释:
物理上,对一个无阻尼或者有阻尼自由振动系统,系统的能量随着时间不应该增加,有阻尼情况还应该减小。因此,一个数值方法的计算结果也不应该放大初始能量,如果经过若干步的数值计算以后,计算结果远比初始条件大,那就是数值算法本身计算是不稳定