编制中心差分法程序报告(结构工程研究生作业)

中心差分法计算单自由度体系的动力反应
北京工业大学结构工程
组员:胡建华 S201204111
马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言
时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。

时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。

由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。

在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。

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

中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为： t
u u u i i ∆-=
-+•
21
1 (a)
2
1
12t
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 时刻的运动方程：
0221
12
11=+∆-+∆+--+-+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 ∆-=
-•
21
10 （h ） `2
1
0102t
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 确定
即 )(1
000ku u c m
u --=•
•
• （k ）
这样就可以根据初始条件00,•
u u 和初始荷载0P ，就可以根据上式确定1-u 的值。

2中心差分法程序编制原理：
（1） 基本数据准备和初始条件计算
)(1
000ku 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 ∆-=
-+•
21
1
2
1
12t
u u u u i i i ∆+-=
-+•
• （4）下一步计算用i+1代替i ，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。

以上为中心差分法逐步计算公式，其具有2阶精度，即误差)(02
t ∆∝ε；并且为有条件稳定，稳定条件为： π
n
T t ≤
∆
上式中，n T 为结构的自振周期，对于多自由度结构体系则为结构的最小自振周期。

3算例
对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量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
2
2
==
=
≤
∆n
n
T t ωπ
所以本次计算取t ∆＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算
4 MATLAB 程序清单
[u,v,ac]=[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 ＝？
5中心差分法计算结果稳定性分析
由以上时程图可以得到当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
2
2
==
=
≤
∆n
n
T t ωπ
，由稳定条件知，当t ∆＝0.4时结果应当是稳定的，而且是发散与收敛的临界点，即t ∆＝0.4为临界步长。

从图中可以看出，此时的计算结果虽然没有发散，但是结果有明显的波动，这预示着若再增大步长t ∆，则计算结果会发散。

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

5
101520-0.05
0.05
时间(s)
位移(m )
5
1015
20
-0.2
-0.100.1
0.2时间(s)
速度(m /s )
05
101520
-2
-10
1时间(s)
加速度(m /s 2)
顶层加速度的时程曲线
5
101520
-0.1
-0.0500.05
0.1时间(s)
位移(m )
顶层位移的时程曲线
5
1015
20
-2-1012x 10
-16时间(s)
速度(m /s )
顶层速度的时程曲线0
5
1015
20
-2-101
2时间(s)
加速度(m /s 2)
顶层加速度的时程曲线
10
203040-0.05
0.05
时间(s)
位移(m )
10
2030
40
-0.2
-0.100.1
0.2时间(s)
速度(m /s )
010
203040
-2
-10
1时间(s)
加速度(m /s 2)
顶层加速度的时程曲线
10
203040
-0.1
-0.0500.05
0.1时间(s)
位移(m )
顶层位移的时程曲线
10
2030
40
-2-1012x 10
-16时间(s)
速度(m /s )
顶层速度的时程曲线0
10
2030
40
-2-101
2时间(s)
加速度(m /s 2)
顶层加速度的时程曲线
10203040
-5
5
时间(s)
位移(m )
010
20
3040
-4
-202
时间(s)
速度(m /s )
010
203040
-5
5
26时间(s)
加速度(m /s 2)
顶层加速度的时程曲线
010
203040
-4-202
48
时间(s)
位移(m )
顶层位移的时程曲线0
10
20
30
40
-2-1012
8时间(s)
速度(m /s )
顶层速度的时程曲线050100
-2
-1012
8
0.405
10
203040-4
-202
时间(s)
位移(m )
10
20
30
40
-1-0.500.5
1
时间(s)
速度(m /s )
010
20
3040
-5
5
5时间(s)
加速度(m /s 2)
顶层加速度的时程曲线。