中心差分法计算程序编程.doc

中心差分法计算程序编程

姓名：张泽伟 学号： 电话：

一、中心差分法程序原理说明

1.1 中心差分法思路

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

1.2 中心差分法原理

中心差分法只在相隔t ?一些离散的时间区间内满足运动方程，其基于有限

差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ?=?，

则速度与加速度的中心差分近似为：

t u u u i i ?-=

-+?211 (a) 21

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 时

刻的运动方程：

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 ）求得体系的速度和加速度。

1.3 初始条件转化

假设给定的初始条件为

),0(),

0(00?

?==u u u u （g ）

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

t u u u ?-=

-?2110 （h ） ` 21

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

确定 即 )(1000ku u c m u --=??

? （k ） 这样就可以根据初始条件

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

值。 1.4 中心差分法编程思路

① 基本数据准备和初始条件计算：

)(1000ku u c m u --=??

? 020012???

-?+?-=u t u t u u ② 计算等效刚度和中心差分计算公式中的相关系数：

t c t m k ?+?=22

22t m

k a ?-

= t c t

m b ?-?=22 ③ 根据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 21

12t u u u u i i i ?+-=-+??

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

1.5 中心差分法稳定条件

以上为中心差分法逐步计算公式，其具有2阶精度，即误差

)(02t ?∝ε；并且为有条件稳定，稳定条件为： πn

T t ≤?

二、程序框图

根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形式展示出来：

三、程序清单

%m,k,c分别为质量、刚度、阻尼

%p0,dt,t分别为外荷载幅值、时间步距、总时间

%u0,v0为初始条件初位移和初速度

%u,v,ac分别为位移、速度、加速度反应

ek=等效刚度；p=荷载；ep=等效荷载

%定义矩阵X0=input('请按格式和顺序输入初始矩阵，如X0=[m,k,c,u0,v0,t,P0,dt],

m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4); %分别取出其中的参数：

v0=X0(1,5);t=X0(1,6);P0=X0(1,7);dt=X0(1,8)

t=[0:dt:t]; %将时间分步，采用等时间步长；

[mm,nn]=size(t); %计算t的向量长度，得出步数；

u=zeros(size(t)); %设定存储u的矩阵；

v=zeros(size(t)); %设定存储v的矩阵；

ac=zeros(size(t)); %设定存储ac的矩阵；

u(:,2)=u0; %赋值向量第2项为u0；

v(:,2)=v0; %赋值向量第2项为v0；

ac(:,2)=(P0-c*v(:,2)-k*u(:,2))/m; %求出初始加速度ac0；

u(:,1)=u(:,2)-dt*v(:,2)+((dt)^2)*ac(:,2)/2; %计算初始条件u-1项；

ek=m/(dt^2)+c/(2*dt); %计算等效刚度；

a=k-(2*m)/(dt^2);

b=m/(dt^2)-c/(2*dt); %计算方程系数；

p(:,2)=P0*sin(0); %给出初始荷载条件；

ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1); %计算初始等效荷载；

u(:,3)=ep(:,2)/ek; %计算位移u1=u(:,3)

for i=3:nn %从第二项开始进行中心差分法计算；

p(:,i)=P0*sin(.5*pi*(i-2)*dt); %给出荷载条件，按照简谐荷载计算；

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

%-----------------------得出所需要结果----------------------------------%

u(:,i+1)=ep(:,i)/ek; %计算位移量；

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

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

end

t=t(:,1:end-1);u=u(:,2:end-1);v=v(:,2:end);ac=ac(:,2:end);p=p(:,2:end);ep=ep(:,2:end);

%------------------------绘制位移、速度、加速度时程曲线-----------------------%

%plot(t,u,'b-o'),hold on,plot(t,v,'g--p'),hold on,plot(t,ac,'r:x'),grid on,xlabel('时间(s)'),ylabel('位移(m)速度(m/s)加速度(m/s^2)'),title('顶层u,v,ac的时程曲线');

subplot(3,1,1),plot(t,u,'b-'),grid,xlabel('时间(s)'),ylabel('位移(m)'),

title('位移u的时程曲线');legend('位移u')

subplot(3,1,2),plot(t,v,'k'),grid,xlabel('时间(s)'),ylabel('速度(m/s)'),

title('速度v 的时程曲线');legend('速度v')

subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),

title('加速度ac 的时程曲线');legend('加速度ac')

四、输入数据

本程序采用单自由度体系进行计算，主要已知参数信息如下：

其质量M=9240kg 、刚度K ＝1460KN/m 、阻尼系数m s kN C /41.6?=，对结构施加动力荷载N t P )5.0sin(73000π=，结构周期T=0.05s ，初始位移

m u 05.00=，初始速度s m v /00=，假设结构处于线弹性状态。由中心差分法可知，要使计算结果稳定且不发散，需满足：时间步长s T t n

159.0=≤?π，

本例分别取时间步长为s 1.0、s 15.0、s 17.0、s 2.0分别进行计算，并验证其稳定条件，取总时间为30s 。则：

X0=[9240 1460000 6410 0.05 0 20 73000 0.05]

五、计算结果

当dt=0.1s:

当dt=0.15s时：

当dt=0.17s时：

当dt=0.2s时：

六、结果稳定性分析

由以上时程图可以得到当t?＝0.1，0.15时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；

当t?＝0.17，0.20时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长t?的增加，计算结果发散得越来越快。

由稳定条件知，当t?≤0.159时结果应当是稳定的，而且是发散与收敛的临界点，所以从以上计算结果可以说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。