《中心差分法计算程序编程.docx》由会员分享，可在线阅读，更多相关《中心差分法计算程序编程.docx（7页珍藏版）》请在金锄头文库上搜索。
1、中心差分法计算程序编程合用标准文案中心差分法计算程序编程姓名：张泽伟学号：电话：一、中心差分法程序原理说明1.1中心差分法思路中心差分法的基本思路：是将运方程中的速度向量和加速度向量用位移的某种合来表示，将微分方程的求解化代数方程的求解，并在区内求得每个渺小区的推公式，而求得整个程的反。1.2中心差分法原理中心差分法只在相隔t一些失散的区内足运方程，其基于有限差分代替位移的求（即速度和加速度），若是采用等步，tit，速度与加速度的中心差分近似：uui1ui12t(a)uui12uiui1t2(b)而失散点的运uiu(ti),uiu(ti),uiu(ti)(i0,1,2,3,)由系统运方程：mu
2、(ti)cu(t)ku(t)0(c)将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以获取ti刻的运方程：mui12uiui1cui1ui1kui0t22t（d）在（d）式中，假ui和ui1是已知的，即在ti及ti以前刻的运已知，可以把已知移到方程的右，整理获取：(mc)ui1(k2m)ui(mc)ui1(e)t22tt2t22t优秀文档合用标准文案由式（e）就可以依照ti及ti以前时刻的运动，求得ti1时刻的运动，若是需要可以用式（a）和式（b）求得系统的速度和加速度。1.3初始条件转变假设给定的初始条件为u0u(0),u0u(0),（g）由式（g）确定u1。在零时刻速度和加速度的
3、中心差分公式为：u1u1u02t（h）u12u0u1u0t2（i）u1u0tu0t2u0（j）将式（i）消去u1得：2而零时刻的加速度值u0可以用t0时的运动方程mu0cu0ku00确定u01(cu0ku0)即m（k）这样就可以依照初始条件u0,u0和初始荷载P0，就可以依照上式确定u1的值。1.4中心差分法编程思路基本数据准备和初始条件计算：u01(cu0ku0)mu1u0tu0t2u02计算等效刚度和中心差分计算公式中的相关系数：mck22tt优秀文档合用标准文案ak2mt2mcbt22t依照ti及ti以前时刻的运动，计算ti1时刻的运动：Pauibui1ui1Pkui1ui1u2tui1
4、2uiui1ut2下一步计算用i+1代替i，对于线弹性结构系统，重复第3步，对于非线性结构系统，重复第2步和第3步。1.5中心差分法牢固条件以上为中心差分法渐渐计算公式，其拥有2阶精度，即误差0(t2)；而且为t有条件牢固，牢固条件为：二、程序框图Tn依照中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形式显现出来：编写回调函数调用输入数据计算初始加速度和位移前一时刻值开始循环体（步数控制）计算等效刚度和系数a、b计算位移、速度、加速度并保存优秀文档合用标准文案解析牢固条件绘图三、程序清单%m,k,c分别为质量、刚度、阻尼%p0,dt,t分别为外荷载幅值、时间步距、总时间%u0,v
5、0为初始条件初位移和初速度%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);%设定储藏
6、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/(dt2)+c/(2*dt);%计算等效刚度；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(
7、:,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)/(dt2);%计算加速胸襟；endt=t(:,
8、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的时程曲线)
9、;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，对结构施加动力荷载P73000sin(0.5t)N
10、，结构周期T=0.05s，初始位移u00.05m，初始速度v00m/s，假设结构处于线弹性状态。由中心差分法可知，要使计算结果牢固且不发散，需满足：时间步长tTn0.159s，本例分别取时间步长为0.1s、0.15s、优秀文档合用标准文案当dt=0.15s时：当dt=0.17s时：优秀文档合用标准文案当dt=0.2s时：六、结果牢固性解析由以上时程图可以获适当t0.1，0.15时渐渐计算结果给出的结构运动趋向收敛的，即计算结果是牢固的；当t0.17，0.20时渐渐计算结果给出的结构运动趋向发散的，即结果是不牢固的，且随着步长t的增加，计算结果发散得越来越快。由牢固条件知，当t0.159时结果应当是牢固的，而且是发散与收敛的临界点，所以从以上计算结果可以说了然中心差分法是有条件牢固的并考据了中心差分法的牢固条件。优秀文档
