中心差分法计算程序编程

中心差分法计算程序编程姓名：张泽伟学号：电话： 一、中心差分法程序原理说明1.1 中心差分法思路中心差分法的基本思路： 是将运动方程中的速度向量和加速度向量用位移的某种组合来表示， 将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。1.2 中心差分法原理中心差分法只在相隔t 一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导 （即速度和加速度） ，如果采用等时间步长，t it ，则速度与加速度的中心差分近似为：uui 12ui 1t(a)uui12uit 2ui 1(b)而离散时间点的运动为uiu(t i ), u iu(ti ), uiu(t i )(i0,1,2,3,)由体系运动方程为：mu (t i )cu(t)ku(t )0(c)将速度和加速度的差分近似公式（a）和式（ b）代入式（ c）可以得到ti 时刻的运动方程：m ui 12uit 2ui 1c ui 12ui 1tkui0（ d）在（d）式中，假设ui 和 ui1 是已知的，即在ti 及t i 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：(mc)u(k2m )u( mc)ut 22ti 1t 2it 22ti 1(e)由式（ e）就可以根据ti 及ti以前时刻的运动，求得t i 1 时刻的运动，如果需要可以用式（ a）和式（ b）求得体系的速度和加速度。1.3 初始条件转化假设给定的初始条件为u0u(0),u 0u(0),（ g）由式（ g）确定 u1 。在零时刻速度和加速度的中心差分公式为：u 0u1u 12t（h）u 0u12u 0u 1t 2（i）u 1u0将式（ i）消去 u1 得：t u 0tu 022（j）而零时刻的加速度值u 0 可以用t 0时的运动方程m u 0c u 0ku00确定u 01 (即mc u 0ku0 )（k）这样就可以根据初始条件u 0 , u0 和初始荷载p0 ，就可以根据上式确定u 1 的值。1.4 中心差分法编程思路 基本数据准备和初始条件计算：u 01 (cu 0mku0 )u 1u 0t u 0tu 022 计算等效刚度和中心差分计算公式中的相关系数：mckt 22t2 makt 2bmct 22t 根据 ti及ti 以前时刻的运动，计算ti 1 时刻的运动：pauibui 1ui 1uuip k1ui 12tuui12uit 2ui 1t)2 下一步计算用 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/(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*(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);%-绘制位移、速度、加速度时程曲线%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/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,xlabel(时间 (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.05m ，初始速度 v00 m/s，假设结构处于线弹性状态。由中心差分法可知，要使计算结果稳定且不发散， 需满足： 时间步长ttn0.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 时逐步计算结果给出