中心差分法计算程序编程.docx
《中心差分法计算程序编程.docx》由会员分享，可在线阅读，更多相关《中心差分法计算程序编程.docx（11页珍藏版）》请在冰豆网上搜索。
中心差分法计算程序编程
中心差分法计算程序编程
姓名：
张泽伟学号：
电话：
一、中心差分法程序原理说明
1.1中心差分法思路
中心差分法的基本思路：
是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
1.2中心差分法原理
中心差分法只在相隔t一些离散的时间区间内满足运动方程，其基于有限
差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，tit，
则速度与加速度的中心差分近似为：
uui1
2
ui1
t
（a）
uui
12ui
t2
ui1
（b）
而离散时间点的运动为
uiu（ti）,ui
u（ti）,ui
u（ti）
（i0,1,2,3,）
由体系运动方程为：
mu（ti）
cu（t）
ku（t）0
（c）
将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到
ti时
刻的运动方程：
mui1
2ui
t2
ui1
cui1
2
ui1
t
kui0
（d）
在（d）式中，假设
ui和ui
1是已知的，即在
ti及ti以前时刻的运动已知，则
可以把已知项移到方程的右边，整理得到：
（mc）u
（k2m）u（mc）u
t22ti1
t2it22ti1
（e）
由式（e）就可以根据
ti及ti
以前时刻的运动，求得
ti1时刻的运动，如果需
要可以用式（a）和式（b）求得体系的速度和加速度。
1.3初始条件转化
假设给定的初始条件为
u0u（0）,
u0u（0）,
（g）
由式（g）确定u
1。
在零时刻速度和加速度的中心差分公式为：
u0u1u1
2t
（h）
u0u1
`
2u0u1
t2
（i）
u1u0
将式（i）消去u1得：
tu0
tu0
2
2
（j）
而零时刻的加速度值
u0可以用t＝0时的运动方程
mu0
cu0
ku00
确定
u01（
即m
cu0
ku0）
（k）
这样就可以根据初始条件
u0,u0和初始荷载
P0，就可以根据上式确定
u1的
值。
1.4中心差分法编程思路
①基本数据准备和初始条件计算：
u0
1（cu0
m
ku0）
u1u0
tu0
tu0
2
2
②计算等效刚度和中心差分计算公式中的相关系数：
mc
kt22t
2m
akt2
bmc
t22t
③根据ti
及ti以前时刻的运动，计算
ti1时刻的运动：
Pauibui1
ui1
uui
Pk
1ui1
2t
uui
12ui
t2
ui1
④下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。
1.5中心差分法稳定条件
以上为中心差分法逐步计算公式，其具有2阶精度，即误差
0（；并且为
有条件稳定，稳定条件为：
tTn
二、程序框图
根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形式展示出来：
编辑回调函数
调用输入数据
三、程序清单
%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））;
ac=zeros（size（t））;
%设定存储v的矩阵；
%设定存储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;
ek=m/（dt^2）+c/（2*dt）;
%计算初始条件u-1项；
%计算等效刚度；
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）
fori=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'）,holdon,plot（t,v,'g--p'）,holdon,plot（t,ac,'r:
x'）,gridon,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、阻尼系数C
6.41kN
s/m，对结
构施加动力荷载P
73000sin（0.5
t）N
，结构周期T=0.05s，初始位移u0
0.05m，
初始速度v0
0m/
s，假设结构处于线弹性状态。
由中心差分法可知，要使计算结果
稳定且不发散，需满足：
时间步长
tTn
0.159s，本例分别取时间步长为
0.1s、0.15s、
0.17s、0.2s分别进行计算，并验证其稳定条件，取总时间为30s。
则：
X0=[9240146000064100.05020730000.05]
五、计算结果
当dt=0.1s:
当dt=0.15s时：
当dt=0.17s时：
当dt=0.2s时：
六、结果稳定性分析
由以上时程图可以得到当t＝0.1，0.15时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；
当t＝0.17，0.20时逐步计算结果给出的结构运动趋向发散的，即结果是
不稳定的，且随着步长t的增加，计算结果发散得越来越快。
由稳定条件知，当t0.159时结果应当是稳定的，而且是发散与收敛的临界点，所以从以上计算结果可以说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。
升级会员 