中心差分法在单自由度中的应用

中心差分法求解单自由度体系的自由振动问题

前言

时域逐步积分法是根据运动方程，引进某些假设，建立由t 时刻状态向量i u 、i u •、i u •

•到t ＋t ∆时刻的状态向量1+i u 、1+•i u 、1+••i u 的递推关系，从而从t ＝0时刻的初始状态向量0u 、0•u 、0•

•u 出发，逐步求出各时刻的状态向量，由于引进的假设条件不同，可以有各种不同的方法，下面主要介绍一种时域逐步积分方法－中心差分法。

中心差分法(central difference method)原理[1]

中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。

中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为：

t

u u u i i ∆-=-+•211 (a) 2

112t u u u u i i i ∆+-=-+•• (b) 而离散时间点的运动为

)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……)

由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c)

将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程： 02211211=+∆-+∆+--+-+i i i i i i ku t u u c t

u u u m （d ） 在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：

12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c t

m u t m k u t c t m (e)

由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。

假设给定的初始条件为

),0(),

0(00••==u u u u （g ）

由式（g ）确定1-u 。在零时刻速度和加速度的中心差分公式为：

t u u u ∆-=

-•2110 （h ） `210102t

u u u u ∆+-=-•• （i ） 将式（i ）消去1u 得：020012

•••-∆+∆-=u t u t u u （j ） 而零时刻的加速度值0••u 可以用t ＝0时的运动方程 0000=++•

••ku u c u m 确定

即 )(1000ku u c m u --=••• （k ） 这样就可以根据初始条件00,•

u u 和初始荷载0P ，就可以根据上式确定1-u 的值。

下面给出采用中心差分法分析时的具体计算步骤：

（1） 基本数据准备和初始条件计算

)(1000ku u c m u --=••• 020012•••

-∆+∆-=u t u t u u （2） 计算等效刚度和中心差分计算公式中的相关系数

t

c t m k ∆+∆=

22 22t

m k a ∆-= t c t m b ∆-∆=22 （3） 根据i t 及i t 以前时刻的运动，计算1+i t 时刻的运动

1---=i i bu au P

k P u i =+1

如果需要，可计算

t u u u i i ∆-=

-+•211 2112t u u u u i i i ∆+-=

-+•• （4）下一步计算用i+1代替i ，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。

以上为中心差分法逐步计算公式，其具有2阶精度，即误差)(02t ∆∝ε；并且为有条件稳定，稳定条件为： πn

T t ≤∆

上式中，n T 为结构的自振周期，对于多自由度结构体系则为结构的最小自振周期。 算例

对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg 、刚度K ＝50KN/m 、阻尼系数C ＝3KNs/m ，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。

采用MATLAB 语言编程，并以单自由度体系为例进行计算,设初位移u0＝0和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件πn

T t ≤∆。

先计算t ∆，由稳定条件n n

T t ωπ2

=≤∆，而52000

50000===M K n ωrad/s, 则4.05

22

===≤∆n n

T t ωπ 所以本次计算取t ∆＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算

MATLAB 程序清单

function [u,v,ac]=centraldifferent(M,C,K,u0,v0,time,dt)

% 本程序采用中心差分法计算结构的动力响应

% 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态；

% ---------%%%%%输入参数%%%%%%%------------

% M------------质量矩阵

% C------------阻尼矩阵

% K------------刚度矩阵

% u0-----------初始位移