中心差分法计算程序编程_第1页
中心差分法计算程序编程_第2页
中心差分法计算程序编程_第3页
中心差分法计算程序编程_第4页
中心差分法计算程序编程_第5页
已阅读5页,还剩3页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

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

2、近似公式(a)和式(b)代入式(c)可以得到时刻的运动方程: (d)在(d)式中,假设和是已知的,即在及以前时刻的运动已知,则可以把已知项移到方程的右边,整理得到: (e)由式(e)就可以根据及以前时刻的运动,求得时刻的运动,如果需要可以用式(a)和式(b)求得体系的速度和加速度。1.3 初始条件转化假设给定的初始条件为 (g)由式(g)确定。在零时刻速度和加速度的中心差分公式为: (h) (i)将式(i)消去得: (j) 而零时刻的加速度值可以用t0时的运动方程 确定即 (k)这样就可以根据初始条件和初始荷载,就可以根据上式确定的值。1.4 中心差分法编程思路 基本数据准备和初始条件计算:

3、计算等效刚度和中心差分计算公式中的相关系数: 根据及以前时刻的运动,计算时刻的运动: 下一步计算用i+1代替i,对于线弹性结构体系,重复第3步,对于非线性结构体系,重复第2步和第3步。1.5 中心差分法稳定条件以上为中心差分法逐步计算公式,其具有2阶精度,即误差;并且为有条件稳定,稳定条件为: 2、 程序框图根据中心差分法的原理,可以得出本程序的主要程序思想,以下面框图的形式展示出来:编辑回调函数调用输入数据计算初始加速度和位移前一时刻值开始循环体(步数控制)计算等效刚度和系数a、b计算位移、速度、加速度并保存分析稳定条件绘图3、 程序清单%m,k,c分别为质量、刚度、阻尼%p0,dt,t分别

4、为外荷载幅值、时间步距、总时间%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); %

5、设定存储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);

6、 %给出初始荷载条件;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(:

7、,i+1)-2*u(:,i)+u(:,i-1)/(dt2); %计算加速度量;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'),hold on,plot(t,v,'g-p'),hold on,plot(t,ac,'r:x'),grid on,xlabel('时间(s)'),ylabel('位移(m)速度(m/s)加速度(m/s2

8、)'),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的时程曲线');legen

9、d('速度v')subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('时间(s)'),ylabel('加速度(m/s2)'),title('加速度ac的时程曲线');legend('加速度ac')4、 输入数据本程序采用单自由度体系进行计算,主要已知参数信息如下:其质量M=9240kg、刚度K1460KN/m、阻尼系数,对结构施加动力荷载,结构周期T=0.05s,初始位移,初始速度,假设结构处于线弹性状态。由中心差分法可知,要使计算结果稳定且不发散,需满足:时间步长,本例分别取时间步长为、分别进行计算,并验证其稳定条件,取总时间为30s。则:X0=9240 1460000 6410 0.05 0 20 73000 0.055、 计算结果当dt=0.1s:当dt=0.15s时:当dt=0.17s时:当dt=0.2s时:6、 结果稳定性分析由以上时程图可以得到当0.1,0.15时逐步计算结果给出的结构运动趋向收敛的,

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

最新文档

评论

0/150

提交评论