于开平-结构动力学第十五讲

( xt t xt ) (1 ) xt (1 )txt t 2
(4) 将计算结果分别赋值给前一时刻，用于下一次的递推计算
xt xt t , xt =xt t , xt xt t
(5) 返回步骤(1)计算下时刻的位移
2.3 纽马克方法(Newmark method)
Mxt t Cxt t Kxt t Qt t
三个方程可以确定三个待定量。由第一个方程可以解出待求时刻加速度表示
xt t =
1 1 1 ( x x ) x ( 1) xt t t t t 2 t t 2
将这个加速度表示代入第二个方程，可以解出待求时刻速度表示
确定步长后，计算系数 a0 (t )2 , a1 2t , a2 (t )2 (4) 计算 xt x0
1 1 x0 x0 2a1 2a0
(5) 形成等效质量阵并对其进行三角分解 M a0 M a1C , • 其中L为三角阵，D为对角阵；
M LDLT
• 等效质量阵一般都是对称的，对称的实矩阵，都可以做这种分解；
中心差分法的几点补充说明 1) 若系统质量矩阵为对角阵，阻尼矩阵也仅建模为与质量阵成比例的形式， 这样等效质量矩阵就是对角阵，其逆矩阵就是对角线元素求倒数，求逆后直
1 接乘上等效力向量，xt t M Qt
就不用三角分解的解代数方程组了，对自由度数多的问题，可大大降低计算 量。这种不需要解代数方程组就可以直接递推求解的方法，称为显示算法。 2) 算法流程中，对时间步长选取的限制，是为了保证数值计算稳定不发散。 为保证数值计算稳定，而限制时间步长大小的算法，称为有条件稳定算法。 ∆������������������ 称为临界步长，根据算法性能分析理论临界步长∆������������������ =
(4) 形成等效刚度矩阵
K K a0 M a1C
(5) 对等效刚度矩阵三角分解
K LDLT
2.3 纽马克方法(Newmark method)
纽马克法的解题步骤——算法流程 2. 对每一时间步进行循环
(1) 计算������ + ∆������ 时刻的等效载荷
Qt t Qt M (a6 xt a2 xt a3 xt ) C (a1 xt a4 xt a5 xt )
������������ ������
= ������
2������ ������ ������
= ������ ，
������
2
������������ 为系统的最小周期，������������ 为系统的最高频率。显然对大型结构，中心差分方 法的时间步长要取得很小。 时间步长越小计算精度越高，但过小，相同计算时间区间内 0, ������ ，需计算的 步数������Τ∆������就大大增加，计算量也跟着增加。因此，时间步长的选择首先要满 足稳定性要求，然后要在精度和计算量之间平衡。
(4) 将计算结果和������时刻的位移分别赋值给前一时刻，用于下一次的递推计算
xt t xt , xt xt t
(5) 返回步骤(1)计算下时刻的位移 根据这个算法流程，同学们可自行用计算机语言来实现。
2.2 中心差分法 (The central difference method)
2.3 纽马克方法(Newmark method)
对待求的下一时刻的位移、速度和加速度在当前时刻������进行泰勒展开
1 1 x (t t ) x (t ) tx (t ) t 2 x (t ) t 3 x (t ) O(t 4 ) 2 6 1 x (t t ) x (t ) tx (t ) t 2 x (t ) O(t 3 ) 2 x (t t ) x (t ) x (tn ) O(t ) x(t t ) x(t ) tx (t ) O(t 2 ) t
将三阶导数表达式代入位移和速度泰勒展开式 t 2 t 2 x (t t ) x (t ) tx (t ) x (t ) x (t t ) O(t 4 ) 3 6
x (t t ) x (t )
t t x (t ) x (t t ) O(t 3 ) 2 2
• 如果不仅对称，同时还正定，可直接做第一讲里介绍的Cholesky分解，
分解成LLT，三角分解后求逆比直接求逆效率要高得多； • 注意第二步零时刻加速度的计算对质量矩阵也要进行这种分解。
2.2 中心差分法 (The central difference method)
中心差分法的解题步骤——算法流程 2. 对每一时间步进行循环 (1) 计算时刻������的等效荷载
x(t ) lim
x(t t ) x(t ) xt t xt t 0 t t
1 x x x ຫໍສະໝຸດ Baidux xt t t t t t t 2 t t
x(t ) lim
x(t ) x(t t ) xt xt t t 0 t t
Mxt + Cxt + Kxt = Qt
2.2 中心差分法 (The central difference method)
将上述中心时刻 t 的速度、加速度表示代入并整理得
1 2 1 K M x M C t xt t 2 2 ( t ) ( t ) 2 t 1 1 2 , a , a 将其中各矩阵前的系数，分别简记为 a0 1 2 (t )2 2t (t ) 2 1 1 M C xt t Qt 2 ( t ) 2 t
直接略去高阶项，并用近似值代替精确值，即
x(t t ) xt t , x(t ) xt , x(t t ) xt t , x(t ) xt , x(t t ) xt t , x(t ) xt
略去高阶项带来的误差，试图用加速度的加权来调节，权系数和保持不变，算法 位移和速度公式为
结 构 动 力 学
第五章 结构动力学中常用的数值算法
（第十五讲）
主讲教师：于开平 哈尔滨工业大学航天学院
2 结构动力响应的数值解法——典型直接积分法
2.1 引言
Mx Cx Kx Q (t ) x (0) x0 x (0) v0
i) 非比例阻尼 ii) 非线性情况F(x,v) iii) 有冲击作用
利用这两种速度表示的平均来确定当前时刻的速度，
1 ( xt t xt t ). 2t
用以当前时刻������为中心的前后时刻位移的差分来计算速度，这也是中心差 分名字的由来。不同时刻的函数值做差，称为函数对时间的差分运算。 加速度用速度的差商表示，每一个时刻的速度再用位移的差商表示
纽马克方法的几点补充说明 1) 常用算法中，为保证计算精度， ������ = 1/2 ，因此也常被称作纽马克 −������ 方法。 为保证计算是无条件稳定的，一般取������ ≥ 1/4。其中������ = 1/4的算法，也被称 为平均加速度方法，相当于在计算的时间区间内假设加速度是不变的。 2) 需要进行一次解代数方程组的运算，不能直接递推求解，无条件稳定的纽 马克−������方法属于隐式算法。
2.4 结构动响应数值算法性能评价指标分析
2.3 纽马克方法(Newmark method)
1 xt t xt txt [( ) xt xt t ]t 2 2 xt t xt [(1 ) xt xt t ]t
其中������, ������ 为加权常数，然后假设������ + ∆������ 时刻的近似值满足运动方程
1) 这三种典型情况，模态叠加不适用，数值积分。 2) 模态分解后，因载荷形式复杂，无法给出解析解，需要数值积分。 这类方法，不需要先进行模态变换，可直接进行积分，因此相比 于模态叠加法，被称为直接积分类方法。
2.2 中心差分法 (The central difference method)
该方法用利用数学上的差商代替导数的思想
Qt Qt ( K a2 M ) xt (a0 M a1C ) xt t
(2) 求解������ + ∆������ 时刻的位移
( LDLT ) xt t Qt
(3) 如需要计算时刻������的速度和加速度值，可通过它们的位移差分表示来计算
xt a1 ( xt t xt t ), xt a0 ( xt t 2 xt xt t )
xt xt t xt 1 xt t xt xt xt t t t t t 1 时间变量在下 ( x 2 x x ) t t t t t 2 标的，表示该 (t )
时刻的近似值
假设上述中心时刻t的速度、加速度连同位移满足振动方程
(2) 确定初始值
x 0 , x 0 , x0
x0 = M 1 (Q (0) Cx0 - Kx0 )
1 1 2
(3) 选择时步长∆������ , 使它满足∆������ < ∆������������������ = ������������ /������(������������ 为系统的最小周期)
(2) 求解������ + ∆������ 时刻的位移
( LDLT ) xt t Qt t
(3) 计算������ + ∆������ 时刻的加速度和速度 1 1 1 xt t = ( x x ) x ( 1) xt t t 2 t t t t 2
xt t
则运动方程可简写成 Mxt t Qt
其中 M a0 M a1C
Qt Qt ( K a2 M ) xt (a0 M a1C ) xt t
x0 a1 ( xt xt ), x0 a0 ( xt 2 x0 xt ) xt =(x0 xt ) / a1
2.3 纽马克方法(Newmark method)
纽马克法的解题步骤——算法流程 1. 初始值计算 (1) 形成系统矩阵K，M和C
(2) 确定初始值 x0 , x0 , x0 ，其中初始加速度的计算与中心差分法一样
(3) 选择时间步长∆t，算法参数γ、β，并计算积分常数
a0
1 1 a 1 a t 2 ， 1 t , 3 2 t a4 1 a5 ( 2) a t (1 ) a t 2 ， , 6 ， 7
xt x0 1 1 x0 x0 2a1 2a0
式中零时刻加速度可通过观测补充给出，也可通过零时刻的动力学方程给出
x0 = M 1 (Q (0) Cx0 - Kx0 )
2.2 中心差分法 (The central difference method)
中心差分法的解题步骤——算法流程 1. 初始值计算 (1) 形成刚度矩阵 K, 质量矩阵 M 和阻尼矩阵C。
xt t
( xt t xt ) (1 ) xt (1 )txt t 2
K K a0 M a1C
将它们同时代入第三个方程，只剩下待求时刻的位移，整理得 Kxt t Qt t
Qt t Qt M (a6 xt a2 xt a3 xt ) C (a1 xt a4 xt a5 xt )