《中心差分法计算程序编程》由会员分享，可在线阅读，更多相关《中心差分法计算程序编程（9页珍藏版）》请在装配图网上搜索。
1、中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明1.1 中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的 某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在 时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。1.2 中心差分法原理中心差分法只在相隔t些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导(即速度和加速度)，如果采用等时间步长， ti t 则速度与加速度的中心差分近似为：uuu i 1 i 1(i 0,1,2,(a)(b)2 tu 2u uUi It2而离散时间点的运动为u
2、 u(t ),u u(t ),u u(t )ii ii ii由体系运动方程为：叫)Cu(t) ku(t) 0(c)将速度和加速度的差分近似公式(a)和式(b)代入式(c)可以得到j时 刻的运动方程：u 2u u u um ic-4z ku 0t22 t i(d)在(d)式中，假设ui和役1是已知的，即在-及-以前时刻的运动已知，则 可以把已知项移到方程的右边，整理得到：m c2m m c()u (k )u ()ut2 2 t i 1t2 it2 2 t i 1(e)由式（e）就可以根据-及S以前时刻的运动，求得Ei时刻的运动，如果需 要可以用式（a）和式（b）求得体系的速度和加速度。1.3 初
3、始条件转化假设给定的初始条件为u u（0）,0u0 u（0）, （g）由式（g）确定u 1。在零时刻速度和加速度的中心差分公式为：u0u0uui12 tu 2u uh)i)uu将式（i）消去ui得：1 0tu0t2u2 0j)而零时刻的加速度值确定0可以用t=0时的运动方程m u 0 cu0ku0u0丄(cum 0ku )0k)10 1t2这样就可以根据初始条件uo，u0和初始荷载P0,就可以根据上式确定u 1的 值。1.4 中心差分法编程思路 基本数据准备和初始条件计算：u07（cu0 ku0）t2u u tuu1 0 0 2 0 计算等效刚度和中心差分计算公式中的相关系数：m ckt2 2
4、 t2ma kt2mcbt2 2 t 根据-及j以前时刻的运动，计算影时刻的运动:P au bui i 1uPki 1uuU -4 I2 tu 2u uUi4i It2 下一步计算用i+1代替i对于线弹性结构体系，重复第3步，对于非线性结 构体系，重复第2 步和第 3 步。1.5 中心差分法稳定条件以上为中心差分法逐步计算公式，其具有2阶精度，即误差（ t2）;并且为T t n 有条件稳定，稳定条件为：二、程序框图根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形 式展示出来：三、程序清单%m,k,c分别为质量、刚度、阻尼%pO,d t,分别为外荷载幅值、时间步距、总时间%u0
5、,v0 为初始条件初位移和初速度%u,v,ac分别为位移、速度、加速度反应 ek二等效刚度；p=荷载；ep=等效荷载%定义矩阵XO=input请按格式和顺序输入初始矩阵，如XO=m,k,c,uO,vO, t,PO,dt, m=XO(l,l);k=XO(l,2);c=XO(l,3);uO=XO(l,%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); u=zeros(size(t); v=zeros(size(t); ac=zeros(size(t);u(:,2)=u0;v(:,2)=v0;
6、 ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m; u(:,l)=u(:,2)-d t*v(:,2) + (d t)22 *ac(:,2)/ ek=m/(d t2)+c/(2*d t);a=k-(2*m)/(d t2);b=m/(d t2)-c/(2*d t);p(:,2)=P0*sin(0);ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1); u(:,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(
7、:,i+1)=ep(:,i)/ek;v:,i)=(u(:,i+1)-u(:,i-1)/(2*dt);% 将时间分步，采用等时间步长；%计算t的向量长度，得出步数；% 设定存储 u 的矩阵；% 设定存储 v 的矩阵；% 设定存储 ac 的矩阵；% 赋值向量第 2 项为 u0；% 赋值向量第 2 项为 v0；% 求出初始加速度 ac0；%计算初始条件u-1项；%计算等效刚度；% 计算方程系数；% 给出初始荷载条件；% 计算初始等效荷载；% 计算位移 u1=u(:,3)% 从第二项开始进行中心差分法计% 给出荷载条件，按照简谐荷载计算% 计算等效荷载；%计算位移量；%计算速度量；ac(:,i) =
8、(u(:,i+l)-2*u(:,i)+u(:,i-l)/(d2); %计算加速度量;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，), hn,dlot(t,v,，gpOhpld)t(t,ac,，r:x，Ongriabe时间(s)，),ylabel(，位移(m)速度(m/s加速度(m/s2), ti顶层u,v,a的时程曲线)；subplot(3,l,l),plot(t,u,b-),gi时间Xiab
9、ei yiab位移(m),title位移u的时程曲线)；legen位移u)subplo t(3,l,2),plo t(t ,v,k),g时间(SabeiylabE度(m/s), titl速度v的时程曲线);legen速度v)subplo t(3,1,3),plo t(t ,ac,r),g时间(SabeiyiabeM 度(m/s2), title加速度ac的时程曲线)；legen加速度ac)四、输入数据本程序采用单自由度体系进行计算，主要已知参数信息如下：其质量M=9240kg、刚度K=1460KN/m、阻尼系数C 6.41kN s/m ,对结构施加动力荷载P 73000 sin(.5 t)N，
10、结构周期T=0.05s初始位移U 血,初始速度。0m/s，假设结构处于线弹性状态。由中心差分法可知，要使计算结果 T稳定且不发散，需满足：时间步长t f 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:(需至7戰貝时间闽当dt二0.15出寸:I! 加速度日Ui1当dt=0.17j时寸:1o加速度日匚1T1051015 时闾筒2025：30位移U1 J位移u的时程曲线5102015 时闾伐） 速威的时程曲线X 1严2530速度汀- 05102015 时间 抑速度就的时程曲线25：30当 dt=0.2S时:-位移U11速度乍11速度”的时程曲线.加速度北11六、结果稳定性分析由以上时程图可以得到当七=01 0.15时逐步计算结果给出的结构运动趋 向收敛的，即计算结果是稳定的；当t = 0.17 0.20时逐步计算结果给出的结构运动趋向发散的，即结果是不 稳定的，且随着步长七的增加，计算结果发散得越来越快。由稳定条件知，当 t 0. 159时结果应当是稳定的，而且是发散与收敛的临 界点，所以从以上计算结果可以说明了中心差分法是有条件稳定的并验证了中心 差分法的稳定条件。
中心差分法计算程序编程
