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

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

2、而离散时间点的运动为? ? ? ui u(ti),ui u(ti ),ui?u(ti )( i0,1,2,3,)? ?由体系运动方程为： mu（ti ） cu（t） ku（t）(c)将速度和加速度的差分近似公式（ a）和式b）代入式（ c）可以得到 ti 时刻的运动方程：d）mui1 2u2i ui1 cui1 ui 1 kui 0t22 t i在（d）式中，假设 ui和ui 1是已知的，即在 ti 及ti以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：( m2c )ui 1 (k 2m2)ui ( m2c )ui 1t2 2 tt 2t 2 2 t (e)由式（ e）就可以根据

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

4、心差分法编程思路 基本数据准备和初始条件计算：?u0( cu0 ku0) m? 2 ? t ?u 1u0tu0u02 计算等效刚度和中心差分计算公式中的相关系数：mct 2 2 t2ma k 2t2mct 2 2 t 根据ti及ti以前时刻的运动，计算 ti 1时刻的运动：Paui bui 1ui 1 P kui 1 ui 12t?uui 1 2ui ui 1 下一步计算用 i+1 代替 i，对于线弹性结构体系，重复第 3 步，对于非线性结 构体系，重复第 2步和第 3步。1.5 中心差分法稳定条件2以上为中心差分法逐步计算公式，其具有 2 阶精度，即误差 0( t2)；并且为有条件稳定，稳定

5、条件为：二、程序框图根据中心差分法的原理， 可以得出本程序的主要程序思想， 以下面框图的形式展示出来：开始循环体(步数控制)计算等效刚度和系数 a、 b计算位移、速度、加速度并保存三、程序清单%m,k,c 分别为质量、刚度、阻尼%p0,dt,t 分别为外荷载幅值、时间步距、总时间 %u0,v0 为初始条件初位移和初速度 %u,v,ac分别为位移、速度、加速度反应 ek=等效刚度； p=荷载； ep=等效荷载%定义矩阵 X0=input(' 请按格式和顺序输入初始矩阵，如m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4); v0=X0(1,5);t=X0(1

6、,6);P0=X0(1,7);dt=X0(1,8) t=0:dt:t;mm,nn=size(t);u=zeros(size(t);v=zeros(size(t);ac=zeros(size(t);u(:,2)=u0;v(:,2)=v0;ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m; u(:,1)=u(:,2)-dt*v(:,2)+(dt)2)*ac(:,2)/ 2; 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

7、(:,3)=ep(:,2)/ek;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); endX0=m,k,c,u0,v0,t,P0,dt, %分别取出其中的参数：%将时间分步，采用等时间步长； %计算 t 的向量长度，得出步数； %设定存储 u 的矩阵； %设定存储 v 的矩阵； %设

8、定存储 ac 的矩阵； %赋值向量第 2 项为 u0； %赋值向量第 2 项为 v0； %求出初始加速度 ac0； %计算初始条件 u-1 项； %计算等效刚度；%计算方程系数； %给出初始荷载条件； %计算初始等效荷载； %计算位移 u1=u(:,3) %从第二项开始进行中心差分法计%给出荷载条件，按照简谐荷载计算； %计算等效荷载；%计算位移量； %计算速度量； %计算加速度量；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);% 绘制位移、速度、加速度时程曲线 %pl

9、ot(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

10、(' 位移 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)'), ti

11、tle(' 加速度 ac 的时程曲线 ');legend(' 加速度 ac')四、输入数据本程序采用单自由度体系进行计算，主要已知参数信息如下：其质量 M=9240kg、刚度 K1460KN/m、阻尼系数 C 6.41kN ? s/ m ，对结构 施加动力荷载 P 73000 sin(0.5 t)N ，结构周期 T=0.05s，初始位移 u0 0.05m ， 初始速度 v0 0m/ s，假设结构处于线弹性状态。 由中心差分法可知，要使计算结果 稳定且不发散， 需满足： 时间步长 t Tn 0.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 时：六、结果稳定性分析由以上时程图可以得到当 t0.1，0.15 时