《中心差分法计算程序编程》由会员分享，可在线阅读，更多相关《中心差分法计算程序编程（9页珍藏版）》请在金锄头文库上搜索。
1、中心差分法计算程序编程姓名：张泽伟学号：电话:一、中心差分法程序原理说明1.1中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。1.2中心差分法原理中心差分法只在相隔厲一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，弋一讥则速度与加速度的中心差分近似为:Ui1Upu=2加(a)Ui1-2UiUiJ(b)而离散时间点的运动为*wi=0,1,2,3,Ui=u(ti),Ui=u
2、(ti),Ui=u(ti)由体系运动方程为：mU（ti）cu（t）ku（t）=0（c）将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到ti时刻的运动方程：Ui1-2554组2Ui1-U.2t(d)在（d）式中，假设Ui和U2是已知的，即在ti及ti以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：)Ui_1(e)(号丄)Ui1(k-卑)Ui-(巴t22t1tt2由式（e）就可以根据1及1以前时刻的运动，求得ti1时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。1.3初始条件转化假设给定的初始条件为Uo二u（0）,UO=U（0）,（g）由式（g）确定
3、u。在零时刻速度和加速度的中心差分公式为：Uo2-t(h)将式（i）消去Ul得:UOU二UoUi-2uou1-：tuO丄UO而零时刻的加速度值U0可以用t=0时的运动方程(i)(j)e*muocuoku0=0确定1UO(-CUo-kUo)m(k)这样就可以根据初始条件Uo,UO和初始荷载Po，就可以根据上式确定的值。 1.4中心差分法编程思路基本数据准备和初始条件计算:1Uo(-CUo-kuO)m t2U4二U0-U0uo2计算等效刚度和中心差分计算公式中的相关系数: 2m:t2c2.：t根据1及ti以前时刻的运动，计算ti1时刻的运动:P二-ambUijiUiIfUi2.：tu二Uii2ui
4、-Uij.：t下一步计算用i+1代替i，对于线弹性结构体系，重复第结构体系，重复第2步和第3步。1.5中心差分法稳定条件3步，对于非线性以上为中心差分法逐步计算公式，其具有2阶精度，即误差MTn有条件稳定，稳定条件为：二二、程序框图根据中心差分法的原理，可以得出本程序的主要程序思想,式展示出来：编辑回调函数I)调用输入数据rs计算初始加速度和位移前一时刻值开始循环体(步数控制)rf勺计算等效刚度和系数a、bL_J计算位移、速度、加速度并保存2乂(M);并且为以下面框图的形三、程序清单%m,k,c分别为质量、刚度、阻尼%pO,dt,t分别为外荷载幅值、时间步距、总时间%uO,vO为初始条件初位移
5、和初速度%u,v,ac分别为位移、速度、加速度反应ek=等效刚度；p=荷载；6卩=等效荷载%t义矩阵XO=input(请按格式和顺序输入初始矩阵，如XO=m,k,c,uO,vO,t,PO,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=O: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)
6、-k*u(:,2)/m;%u(:,1)=u(:,2)-dt*v(:,2)+(dt)A2)*ac(:,2)/2;%ek=m/(dtA2)+c/(2*dt);%a=k-(2*m)/(dtA2);b=m/(dtA2)-c/(2*dt);%p(:,2)=P0*sin(0);%ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1);%u(:,3)=ep(:,2)/ek;%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
7、(:,i)=(u(:,i+1)-u(:,i-1)/(2*dt);%ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1)/(dtA2);%分别取出其中的参数：将时间分步，采用等时间步长；计算t的向量长度，得出步数；设定存储u的矩阵；设定存储v的矩阵；设定存储ac的矩阵；赋值向量第2项为u0；赋值向量第2项为v0；求出初始加速度ac0；计算初始条件u-1项;计算等效刚度；计算方程系数；给出初始荷载条件；计算初始等效荷载；计算位移u仁u(:,3)从第二项开始进行中心差分法计给出荷载条件，按照简谐荷计算等效荷载；%计算位移量；计算速度量；计算加速度量；endt=t(:,1:end-1
8、);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/sA2),title(顶层u,v,ac时间(s),ylabel(位移(m),位移u)时间(s),ylabel(速度(m/s),速度v)时间(s),ylabel(加速度的时程曲线);subplot(3,1,1),plo
9、t(t,u,b-),grid,xlabel(title(位移u的时程曲线);legend(subplot(3,1,2),plot(t,v,k),grid,xlabel(title(速度v的时程曲线);legend(m/sA2),title(subplot(3,1,3),plot(t,ac,r),grid,xlabel(加速度ac的时程曲线）;legend（加速度ac）四、输入数据本程序采用单自由度体系进行计算，主要已知参数信息如下:其质量M=9240kg刚度K=1460KN/m阻尼系数C=6.41kN*s/m，对结构施加动力荷载P=73000sin（05：t）N，结构周期T=0.05s，初始位
10、移Uo05m,初始速度V。=0m/s，假设结构处于线弹性状态。由中心差分法可知，要使计算结果稳定且不发散，需满足：时间步长T=0.159s，本例分别取时间步长为0.1s、0.15s、0.17s、0.2s分别进行计算，并验证其稳定条件，取总时间为30s。则：X0=9240146000064100.05020730000.05五、计算结果当dt=0.1s:时间速度的时程曲线hl:LL加速度敢1II51015时间202630时间加速度貳的时程曲线100当dt=0.15s时:时间速度的时程曲线052530101520时间加速度兀的时程曲线P-i1加速度日匚LIII5001015202530时间（旬60
11、4020当dt=0.17s时:50-位移u的时程曲线Iriir位移uiiI13CI101620时间速度甲的时程曲线L速度訂IIiri255101520时间加速度能的时程曲线1iIi加速度ac|iiTT101520时间当dt=0.2s时:8I1i位移uLi1X105位移u的时程曲线101520时间（旬速度/的时程曲线25300-加速度礼1iiix10C751015时间20263011!Dn速度常I1ii5101620时间加速度机的时程曲线2530六、结果稳定性分析由以上时程图可以得到当t=0.1,0.15时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；当氏=0.17,0.20时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长人t的增加，计算结果发散得越来越快。由稳定条件知，当人t0.159时结果应当是稳定的，而且是发散与收敛的临界点，所以从以上计算结果可以说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。欢迎您的下载，资料仅供参考!致力为企业和个人提供合同协议，策划案计划书，学习资料等等打造全网一站式需求
