中心差分法

位 移 (m)
2 0 0 x 10
-3
0.5
1
1.5
速 度 (m/s)
5 0 -5
2.5 3 时 间 (s) 速 度 v的 时 程 曲 线
2
3.5
4
4.5
5
速度v
0
0.5
1
1.5
加 速 度 (m/s 2)
0.5 0 -0.5 0 0.5 1 1.5
2.5 3 3.5 时 间 (s) 加 速 度 ac 的 时 程 曲 线
%线弹性条件下，可计算单自由度或多自由度振动问题 %-----------输入参数-------------------。 %M 质量矩阵，C 阻尼矩阵，K 刚度矩阵，u0 初始时刻位移矩阵，v0 初始时刻速度矩阵。time 求解的时 间区间，dt 时间步长。 M=input('输入质量矩阵 M ：'); C=input('输入 C 阻尼矩阵 C：'); K=input('输入刚度矩阵 K：'); u0=input('输入初始位移矩阵 u0：'); v0=input('输入初始速度矩阵 v0:'); time=input('输入求解时间区间'); dt=input('输入求解步长 dt:'); m=input('输入求解自由度数 m：'); P0=input('输入荷载幅值 P0：'); %u 输出位移矩阵，v 输出速度矩阵，ac 输出加速度矩阵。 %----------求解值---------------------。 %设定存储矩阵 t=[0:dt:time]; n=time/dt; u=zeros(m,n+1); v=zeros(m,n+1); ac=zeros(m,n+1); P=zeros(m,n+1); u(:,2)=u0; v(:,2)=v0; %-----计算等效系数--------------------。 Ks=M/dt^2+C/(2*dt); a=K-2*M/dt^2; b=M/dt^2-C/(2*dt); %---------进行计算--------------。 ac(:,2)=P0/M-1/M*(C*v(:,2)+K*u(:,2)); u(:,1)=u(:,2)-dt*v(:,1)+dt^2/2*ac(:,2); u(:,3)=(P0 -a*u(:,2)+b*u(:,1))/Ks; for i=3:n P(:,i)=P0*sin(6*pi*(i-2)*dt);
假设 u i , u i 1 ，已知，把已知项整理到右边，
(
m c m c ) u i 1 F (t ) (k 2m ) u i ( ) ui 1 2 2 2 t 2 t t t t 2
这样就可以根据 t i , t i 1 求出 t i 1 运动参数。 二，中心差分法实现基本过程 （1），计算初始条件：
u(:,i+1)=(P(:,i)-a*u(:,i)-b*u(:,i-1))/Ks; 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);v=v(:,2:end);ac=ac(:,2:end); 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') ;
0 -2 0 x 10
-3
0.5
1
1.5
速 度 (m/s)
5 0 -5
2.5 3 时 间 (s) 速 度 v的 时 程 曲 线
2
3.5
4
4.5
5
速度v
0
0.5
1
1.5
加 速 度 (m/s 2)
0.5 0 -0.5 0 0.5 1 1.5
2.5 3 3.5 时 间 (s) 加 速 度 ac 的 时 程 曲 线
u 0 u (0); . . u 0 u (0);
. u1 u 1 u0 2t ；
.. u1 2 u 0 u 1 ； u0 t 2
. .. t 2 消去 u1 ，得到 u 1 u 0 t u 0 2 u 0 ；
.. . 1 u 0 P 0 / m (c u 0 k u 0) 。 m
0
0.5
1
1.5
加 速 度 (m/s 2)
0.5 0 -0.5 0 0.5 1 1.5
2.5 3 3.5 时 间 (s) 加 速 度 ac 的 时 程 曲 线
2
4
4.5
5
加 速 度 ac
2
2.5 时 间 (s)
3
3.5
4
4.5
5
t =0.1,
2 x 10
-3
位 移 u的 时 程 曲 线 位移u
位 移 (m)
（4），对于线弹性问题，重复步骤（3）。 三，中心差分法计算实例。 振动系统质量为 10000kg,阻尼为 10KN∙s/m，K＝,1500000N/m，外 载荷为 P=5000∙sin6πN。初始位移和速度均为零。 1，稳定性分析： 中心差分法的稳定条件：t
Tn 2 2 m =0.164. k n
2
4
4.5
5
加 速 度 ac
2
2.5 时 间 (s)
3
3.5
4
4.5
5
t =0.05,
5 x 10
-4
位 移 u的 时 程 曲 线 位移u
位 移 (m)
0 -5 0 x 10
-3
0.5
1
1.5
速 度 (m/s)
5 0 -5
2.5 3 时 间 (s) 速 度 v的 时 程 曲 线
2
3.5
4
4.5
5
速度v
. u i 1 ui 1 ui 2t
.. u i 1 2u i u i 1 u t 2
振动方程为：
.. . m u c u ku F (t )
将加速度与速度的表达式带入：
2u i u i 1 u u u m i 1 c i 1 i 1 ku F 2t t 2
中心差分法 一，中心差分法基本原理： 1，将运动方程中的速度向量和加速度向量用位移的某种组合来 表示，将微分方程组的求解问题转化为代数方程组的求解问题，并 在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反 应,。 中心差分法是一种显示积分法，用位移的有限差分代替位移的导数即加速 度和速度。
速度v
0
0.5
1
1.5
加 速 度 (m/s 2)
0.5 0 -0.5 0 0.5 1 1.5
2.5 3 3.5 时 间 (s) 加 速 度 ac 的 时 程 曲 线
2
4
4.5
5
加 速 度 ac
2
2.5 时 间 (s)
3
3.5
4
4.5
5
比较以上结果，时间步长越小，精度越高，当不符合稳定条件后，后 接近不稳定条件大约 1/2 时，结果就十分的粗犷，精度很低。加速 度由于解法自身缺陷问题，在零时刻附近震荡剧烈，不准确。 附录： Matlab 求解代码：
2
3.5
4
4.5
5
速度v
0
0.5
1
1.5
加 速 度 (m/s 2)
0.5 0 -0.5 0 0.5 1 1.5
2.5 3 3.5 时 间 (s) 加 速 度 ac 的 时 程 曲 线
2
4
4.5
5
加 速 度 ac
2
2.5 时 间 (s)
3
3.5
4
4.5
5
t =0.01,
4 x 10
-4
位 移 u的 时 程 曲 线 位移u
（2），计算等效刚度及中心差分公式当中的系数：
m c k ; t 2 2t
ak 2m ; t 2
b
m c 。 t 2 2t
（3），根据 t i , t i 1 ，求出 t i 1 的运动参数
ui 1
F (t ) a u i b u i 来自 1 _ k24
4.5
5
加 速 度 ac
2
2.5 时 间 (s)
3
3.5
4
4.5
5
t =0.2,
0 位 移 u的 时 程 曲 线 位移u
位 移 (m)
-0.005 -0.01 0 x 10
-3
0.5
1
1.5
速 度 (m/s)
5 0 -5
2.5 3 时 间 (s) 速 度 v的 时 程 曲 线
2
3.5
4
4.5
5
则t =0.005,0.01，0.05,0.1，0.2 分别进行计算。 2，计算结果如下：
t =0.005
4 x 10
-4
位 移 u的 时 程 曲 线 位移u
位 移 (m)
2 0 0 x 10
-3
0.5
1
1.5
速 度 (m/s)
5 0 -5
2.5 3 时 间 (s) 速 度 v的 时 程 曲 线