中心差分法计算程序编程

中心差分法计算程序编程 姓名：张泽伟 学号： 电话： 一、中心差分法程序原理说明 1.1 中心差分法思路 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移 的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题， 并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 1.2 中心差分法原理 中心差分法只在相隔 一些离散的时间区间内满足运动方程，其基于有限t 差分代替位移对时间的求导（即速度和加速度） ，如果采用等时间步长， ，则速度与加速度的中心差分近似为：ti (a)t uii21 (b)2 11tiii 而离散时间点的运动为 ( 0,1,2,3,),(),(iiii tutut 由体系运动方程为： (c)0kcmi 将速度和加速度的差分近似公式（a）和式（b）代入式（ c）可以得到 时it 刻的运动方程： （d） 022111 iiiiii kutctu 在（d）式中，假设 和 是已知的，即在 及 以前时刻的运动已知，i1i iti 则可以把已知项移到方程的右边，整理得到： (e)12212 )()()( iii utctmutkutctm 由式（e）就可以根据 及 以前时刻的运动，求得 时刻的运动，如果需iit 1i 要可以用式（a）和式（b ）求得体系的速度和加速度。 1.3 初始条件转化 假设给定的初始条件为 （g）),0( u 由式（g）确定 。在零时刻速度和加速度的中心差分公式为：1u （h）t 210 （i）2 1010tuu 将式（i）消去 得： （j） 1 001t 而零时刻的加速度值 可以用 t0 时的运动方程 0 u 00kucm 确定 即 （k） )(100kucm 这样就可以根据初始条件 和初始荷载 ，就可以根据上式确定 的0, 0P1u 值。 1.4 中心差分法编程思路 基本数据准备和初始条件计算： )(100kucmu0201t 计算等效刚度和中心差分计算公式中的相关系数： tctmk2atctb2 根据 及 以前时刻的运动，计算 时刻的运动：iti 1i1iiuaPki1tuii21211tiii 下一步计算用 i+1 代替 i，对于线弹性结构体系，重复第 3 步，对于非线性 结构体系，重复第 2 步和第 3 步。 1.5 中心差分法稳定条件 以上为中心差分法逐步计算公式，其具有 2 阶精度，即误差 ；并且)(02t 为有条件稳定，稳定条件为： nTt 2、程序框图 根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的 形式展示出来： 编辑回调函数 调用输入数据 计算初始加速度和位移前一时刻值 计算等效刚度和系数 a、 b 开始循环体（步数控制） 3、程序清单 %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); %设定存储 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(:,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); %计算加速度量； 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),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,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) 4、输入数据 本程序采用单自由度体系进行计算，主要已知参数信息如下： 其质量 M=9240kg、刚度 K1460KN/m、阻尼系数 ，对结mskNC/41.6 构施加动力荷载 ，结构周期 T=0.05s，初始位移 ，NtP)5.0sin(73 u05. 初始速度 ，假设结构处于线弹性状态。由中心差分法可知，要使计算结果mv/0 稳定且不发散，需满足：时间步长 ，本例分别取时间步长为 、sTtn19.s1.0 、 、 分别进行计算，并验证其稳定条件，取总时间为 30s。则：s15.07.s2. X0=9240 1460000 6410 0.05 0 20 73000 0.05 5、计算结果 当 dt=0.1s: 当 dt=0.15s 时： 当 dt=0.17s 时： 当 dt=0.2s 时： 6、结果稳定性分析 由以上时程图可以得到当 0.1，0.15 时逐步计算结果给出的结构运动趋t 向