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

1、中心差分法计算程序编程中心差分法计算程序编程姓名：张泽伟 学号： 电话： 一、中心差分法程序原理说明1.1 中心差分法思路中心差分法的基本思路： 是将运动方程中的速度向量和加速度向量用位移的某种组合来表示， 将微分方程组的求解问题转化为代数方程组的求解问题， 并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。1.2 中心差分法原理中心差分法只在相隔 t 一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导 （即速度和加速度） ，如果采用等时间步长， t i t ，则速度与加速度的中心差分近似为：u ui 12ui 1t(a)u ui1 2uit 2ui 1(

2、b)而离散时间点的运动为ui u(t i ), u iu(ti ), uiu(t i )( i 0,1,2,3, )由体系运动方程为：mu (t i )cu(t)ku(t ) 0(c)将速度和加速度的差分近似公式（ a）和式（ b）代入式（ c）可以得到ti 时刻的运动方程：m ui 12uit 2ui 1c ui 12ui 1tkui 0（ d）在（d）式中，假设ui 和 ui1 是已知的，即在ti 及t i 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：( m c )u(k 2m )u ( m c )ut 2 2 t i 1t 2 i t 2 2 t i 1(e)由式（ e）

3、就可以根据ti 及ti以前时刻的运动，求得t i 1 时刻的运动，如果需要可以用式（ a）和式（ b）求得体系的速度和加速度。1.3 初始条件转化假设给定的初始条件为u0 u(0),u 0 u(0),（ g）由式（ g）确定 u1 。在零时刻速度和加速度的中心差分公式为：u 0 u1 u 12 t（h）u 0 u12u 0 u 1t 2（i）u 1 u0将式（ i）消去 u1 得：t u 0t u 022（j）而零时刻的加速度值u 0 可以用 t 0 时的运动方程m u 0c u 0ku0 0确定u 0 1 (即 mc u 0ku0 )（k）这样就可以根据初始条件u 0 , u0 和初始荷载P

4、0 ，就可以根据上式确定u 1 的值。1.4 中心差分法编程思路 基本数据准备和初始条件计算：u 01 ( cu 0mku0 )u 1 u 0t u 0t u 022 计算等效刚度和中心差分计算公式中的相关系数：m ck t 2 2 t2 ma k t 2b m ct 2 2 t 根据 ti及ti 以前时刻的运动，计算ti 1 时刻的运动：P aui bui 1ui 1u uiP k1 ui 12 tu ui1 2uit 2ui 1 下一步计算用 i+1 代替 i，对于线弹性结构体系，重复第 3 步，对于非线性结构体系，重复第 2 步和第 3 步。1.5 中心差分法稳定条件以上为中心差分法逐步

5、计算公式，其具有 2 阶精度，即误差0( ；并且为有条件稳定，稳定条件为：t Tn二、程序框图根据中心差分法的原理， 可以得出本程序的主要程序思想， 以下面框图的形式展示出来：编辑回调函数调用输入数据三、程序清单%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,

6、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)/

7、m;% 求出初始加速度 ac0 ；u(:,1)=u(:,2)-dt*v(:,2)+(dt)2)*ac(:,2)/2;ek=m/(dt2)+c/(2*dt);% 计算初始条件 u-1 项；% 计算等效刚度；a=k-(2*m)/(dt2);b=m/(dt2)-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*

8、(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)/(dt2); % 计算加速度量； endt=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);%-

9、 绘制位移、速度、加速度时程曲线 %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/s2),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,xlab

10、el( 时间 (s),ylabel( 速度 (m/s), title( 速度 v 的时程曲线 );legend( 速度 v)subplot(3,1,3),plot(t,ac,r),grid,xlabel( 时间 (s),ylabel( 加速度 (m/s2), title( 加速度 ac 的时程曲线 );legend( 加速度 ac)四、输入数据本程序采用单自由度体系进行计算，主要已知参数信息如下：其质量 M=9240kg 、刚度 K1460KN/m 、阻尼系数 C6.41kNs/ m ，对结构施加动力荷载 P73000 sin( 0.5t )N，结构周期 T=0.05s ，初始位移 u00.0

11、5m ，初始速度 v00 m/s，假设结构处于线弹性状态。 由中心差分法可知，要使计算结果稳定且不发散， 需满足： 时间步长t Tn0.159s ，本例分别取时间步长为0.1s 、0.15s 、0.17s 、 0.2s分别进行计算，并验证其稳定条件，取总时间为 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 时结果应当是稳定的，而且是发散与收敛的临界点，所以从以上计算结果可以说明了中心差分法是有条件稳定的并验证了中心 差分法的稳定条件。