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

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

前言

时域逐步积分法是根据运动方程，引进某些假设，建立由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-----------初始位移

% v0-----------初始速度

% time---------模拟时间

% dt-----------时间步长

% ---%%%%%%输出值%%%%%%%%------

% u--------------位移

% v--------------速度

% ac-------------加速度

% -------%%%%%%%%中心差分法主要公式及原理%%%%%%%%%%----------- % MX''+CX'+KX=0

% M*(X(t+dt)-2*X(t)+X(t-dt))/(dt^2)+C*(X(t+dt)-X(t-dt))/(2*dt)+K*X(t)=0

% (M/dt^2+C/2*dt)*(X(t+dt))=-(K-2*M/dt^2)*X(t+dt)-(M/dt^2-C/2*dt)*X(t-dt)

%----------------- 等效刚度Ke等效荷载Pe和相关系数a,b-------------------------

% Ke=M/dt^2+C/2*dt

% a=K-2*M/dt^2

% b=M/dt^2-C/2*dt

% Pe=-a*X(t)-b*X(t-dt)

% X(t+dt)=Pe/Ke

% X(t)'=(X(t+dt)-X(t-dt))/(2*dt)

% X(t)''=(X(t+dt)-2*X(t)+X(t-dt))/(dt^2)

% ------------------初始条件---------------------

% X0''=(-C*X0'-K*X0)

% X(-1)=X0-X0'*dt+X0''*(dt^2)/2

% --------@Copyright by zhouhuaping(S201004232)-----

clear all

M=input('输入质量矩阵M ：');

C=input('输入阻尼矩阵C：');

K=input('输入刚度矩阵K：');

u0=input('输入初始位移u0：');

v0=input('输入初始速度v0：');

time=input('输入模拟时间time：');

dt=input('输入时间步长dt ：');

[m,m]=size(K);

n=time/dt; %计算步数

u=zeros(m,floor(n)+1); %设定存储位移矩阵

v=zeros(m,floor(n)+1); %设定存储速度矩阵

ac=zeros(m,floor(n)+1); %设定存储加速度矩阵

P=zeros(m,floor(n)+1); %设定存储荷载矩阵

u(:,2)=u0; %给定初位移

v(:,2)=v0; %给定初速度

Ke=M/(dt^2)+((C)/(2*dt)); %等效刚度Ke及系数a、b a=K-2*M/dt^2;

b=M/dt^2-C/(2*dt);

for i=3:1:floor(n)+1;

t=(i-2)*dt;

ac(:,2)=M\(-K*u(:,2)-C*v(:,2)); %计算初加速度

u(:,1)=u(:,2)-v(:,2)*dt+(ac(:,2)*(dt^2))/2; %计算(0-dt)时刻位移

Pe= -a*u(:,i-1)-b*u(:,i-2); %计算等效荷载Pe

u(:,i)=Ke\Pe; %计算位移

v(:,i)=(u(:,i)-u(:,i-2))/(2*dt); %计算速度

ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2))/(dt^2); %计算加速度

end

%--------%%%%%%%%%%绘制位移、速度、加速度时程曲线%%%%%%%%%%%-------- t=0:dt:time;

subplot(2,2,1),plot(t,u(m,:),'k-'),grid,xlabel('时间(s)'),ylabel('位移(m)'),title('顶层位移的时程曲线');

subplot(2,2,2),plot(t,v(m,:),'r-'),grid,xlabel('时间(s)'),ylabel('速度(m/s)'),title('顶层速度的时程曲线');

subplot(2,2,3),plot(t,ac(m,:),'b-'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),title('顶层加速度的时程曲线');

%----------end

运行centraldifferent.M 文件

输入参数：

K=50000; M=2000; C=3000; u0=0; v0=0;time ＝20s ；dt ＝？

中心差分法计算结果稳定性分析

由以上时程图可以得到当t ?＝0.1, 0.3, 0.4时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；当t ?＝0.41,0.42, 0.45时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长t ?的增加，计算结果发散得越来越快。 由稳定条件n n

T t ωπ2

=≤?，而52000

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

22

===≤?n n

T t ωπ，由稳定条件知，当t ?＝0.4时结果应当是稳定的，而且是发散与收敛的临界点，即t ?＝0.4为临界步长。从图中可以看出，此时的计算结果虽然没有发散，但是结果有明显的波动，这预示着若再增大步长t ?，则计算结果会发散。

以上计算分析，说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。