结构动力学使用中心差分法计算单自由度体系动力反应的MATLAB程序.docx

《结构动力学使用中心差分法计算单自由度体系动力反应的MATLAB程序.docx》由会员分享，可在线阅读，更多相关《结构动力学使用中心差分法计算单自由度体系动力反应的MATLAB程序.docx（8页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上中心差分法计算单自由度体系动力反映的报告前言基于叠加原理的时域积分法与频域积分法一样，都假设结构在在全部反应过程中都是线性的。而时域逐步积分法只是假设结构本构关系在一个微小的时间步距内是线性的，相当于分段直线来逼近实际的曲线。时域逐步积分法是结构动力问题中研究并应用广泛的课题。中心差分法是一种目前发展的一系列结构动力反应分析的时域逐步积分法的一种，时域逐步积分法还包括分段解析法、平均常加速度法、线性加速度法、Newmarket-和Wilson-法等。中心差分法(central difference method)原理 中心差分法的基本思路将运动方程中的速度向量和加速度

2、向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。中心差分法是一种显示的积分法，它基于用有限差分代替位移对时间的求导（即速度和加速度）。如果采用等时间步长，ti=t（t为常数），则速度与加速度的中心差分近似为ui=ui+1+ui-12t (1)ui=ui+1-2ui+ui-1t2 (2)用u表示位移，离散时间点的运动为：ui=uti，ui=uti，ui=uti (i=0,1,2)体系的运动方程为mut+cut+kut=P(t) (3)将速度和加速度的差分近似公式(1)和(2)代入(3)中得出在ti

3、时刻的运动方程，将方程整理得到ui+1由ui和ui-1表示的两步法的运动方程(4)：mt2+c2tui+1=Pi-k-2mt2ui-mt2-c2tui-1 (4)这样就可以根据ti及以前的时刻的运动求得ti+1时刻的运动。中心差分法属于两步法，用两步法计算时存在起步问题，必须要给出相邻的两个时刻的位移值，才能逐步计算。对于地震作用下结构的反应问题和一般的零初始条件下的动力问题，可以用(4)直接计算，因为总可以假设初始的两个时间点（一般取i=0,-1）的位移等于零。但是对应于非零初始条件或零时刻外荷载很大时，需要进行一定的分析，建立两个起步时刻（即i=0,-1）的位移值。假设给定的初始条件为u0

4、=u0u0=u0 (5)根据初始条件来确定u-1。根据中心差分公式u0=u1+u-12tu0=u1-2u0+u-1t2 (6)消去u1得到u-1的公式：u-1=u0-tu0+t22u0 (7)其中零时刻加速度值u0可以由t=0时的运动方程得到即u0=1mP0-cu0-ku0 (8)这样就可以根据初始条件得到u-1，然后再将初始条件应用于公式(4)中，逐步求出不同时刻的运动。中心差分法分析时的具体计算步骤：(1) 基本数据准备与初始条件计算已知：初始位移u0、u0和初始荷载值P0来计算u0和u-1u0=1mP0-cu0-ku0u-1=u0-tu0+t22u0 (2) 计算等效刚度和中心差分法计算

5、公式中的系数k=mt2+c2ta=k-2mt2b=mt2-c2t因此中心差分法计算公式可以表示为：kui+1=Pi-aui-bui-1(3) 根据ti及以前的时刻的运动求得ti+1时刻的运动Pi=Pi-aui-bui-1ui+1=Pik(4)下一步计算中用i+1代替i，对于线弹性体系重复第3步计算步骤，对于非线性弹性体系，重复第2和第3计算步骤。以上的中心差分法逐步计算公式具有2阶精度，即误差O(t2)；并且是有条件稳定的，稳定条件为：tTn式中，Tn为结构的自振周期，对于多自由度体系则为结构的最小自振周期。算例本算例根据结构动力学48页算例3.1数据编写，稳定条件为dt=0.16s对于一个单

6、层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=9240kg、刚度K1460KN/m、阻尼系数C6.41KNs/m，对结构施加动力荷载P=73000sin0.5t假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u00.05m和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件。先计算，由稳定条件，而rad/s，则，所以本次计算取0.2，0.1，0.05分别进行计算。计算结果与分析 1）当dt=0.05s时，可以得到位移u，速度v，加速度ac的时程曲线如下：2）当dt=01s时，可以得到位移

7、u，速度v，加速度ac的时程曲线如下：3）当dt=0.2s时，可以得到如下提示：不满足稳定条件:dt2*sqrt(m/k) %判断时间步长是否满足稳定条件 disp(不满足稳定条件:dt=Tn/pi,请重新输入符合稳定条件的时间步长dt) return elseif 0dt=2*sqrt(m/k) disp(满足稳定条件为：dt=Tn/pi) endt=0:dt:all_time; %将时间分步，采用等时间步长；mm,nn=size(t); %计算t的向量长度，得出步数；u=zeros(size(t); %设定存储u的矩阵；v=zeros(size(t); %设定存储v的矩阵；ac=zeros

8、(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);

9、 %计算初始等效荷载；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);%计算加速度量；

10、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),title(顶层u,v,ac的时程曲线);subplot(3,1,1),plot(t,u,b-),grid,xlabel(时间(s),ylabel(位移(m),tit

11、le(位移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)中心差分法求解单自由度体系的自由振动问题算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg、刚度K50KN/m、阻尼系数C3KNs/m，假设结构处于线

12、弹性状态，用中心差分法计算结构的自由振动反应。采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u00和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件。先计算，由稳定条件，而rad/s,则所以本次计算取0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算MATLAB程序清单function u,v,ac=centraldifferent(M,C,K,u0,v0,time,dt) % 本程序采用中心差分法计算结构的动力响应% 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态；% -%输入参数%-% M-质量

13、矩阵% C-阻尼矩阵% K-刚度矩阵% u0-初始位移% v0-初始速度% time-模拟时间% dt-时间步长% -%输出值%-% u-位移% v-速度% ac-加速度% -%中心差分法主要公式及原理%-% MX+CX+KX=0 % M*(X(t+dt)-2*X(t)+X(t-dt)/(dt2)+C*(X(t+dt)-X(t-dt)/(2*dt)+K*X(t)=0% (M/dt2+C/2*dt)*(X(t+dt)=-(K-2*M/dt2)*X(t+dt)-(M/dt2-C/2*dt)*X(t-dt)%- 等效刚度Ke等效荷载Pe和相关系数a,b-% Ke=M/dt2+C/2*dt% a=K-

14、2*M/dt2% b=M/dt2-C/2*dt% Pe=-a*X(t)-b*X(t-dt)% X(t+dt)=Pe/Ke% X(t)=(X(t+dt)-X(t-dt)/(2*dt)% X(t)=(X(t+dt)-2*X(t)+X(t-dt)/(dt2)% -初始条件-% X0=(-C*X0-K*X0)% X(-1)=X0-X0*dt+X0*(dt2)/2% -Copyright by zhouhuaping(S)-clear allM=input(输入质量矩阵M ：);C=input(输入阻尼矩阵C：);K=input(输入刚度矩阵K：);u0=input(输入初始位移 u0： );v0=in

15、put(输入初始速度 v0： );time=input(输入模拟时间 time：);dt=input(输入时间步长dt ：);m,m=size(K);n=time/dt; %计算步数u=zeros(m,floor(n)+1); %设定存储位移矩阵 v=zeros(m,floor(n)+1); %设定存储速度矩阵 ac=zeros(m,floor(n)+1); %设定存储加速度矩阵 P=zeros(m,floor(n)+1); %设定存储荷载矩阵 u(:,2)=u0; %给定初位移v(:,2)=v0; %给定初速度Ke=M/(dt2)+(C)/(2*dt); %等效刚度Ke及系数a、ba=K-2

16、*M/dt2;b=M/dt2-C/(2*dt);for i=3:1:floor(n)+1; t=(i-2)*dt; ac(:,2)=M(-K*u(:,2)-C*v(:,2); %计算初加速度 u(:,1)=u(:,2)-v(:,2)*dt+(ac(:,2)*(dt2)/2; %计算(0-dt)时刻位移 Pe= -a*u(:,i-1)-b*u(:,i-2); %计算等效荷载Pe u(:,i)=KePe; %计算位移 v(:,i)=(u(:,i)-u(:,i-2)/(2*dt); %计算速度 ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2)/(dt2); %计算加速度end%

17、-%绘制位移、速度、加速度时程曲线%- t=0:dt:time; subplot(2,2,1),plot(t,u(m,:),k-),grid,xlabel(时间(s),ylabel(位移(m),title(顶层位移的时程曲线); subplot(2,2,2),plot(t,v(m,:),r-),grid,xlabel(时间(s),ylabel(速度(m/s),title(顶层速度的时程曲线); subplot(2,2,3),plot(t,ac(m,:),b-),grid,xlabel(时间(s),ylabel(加速度(m/s2),title(顶层加速度的时程曲线); %-end运行centra

18、ldifferent.M文件输入参数：K=50000; M=2000; C=3000; u0=0; v0=0;time20s；dt？中心差分法计算结果稳定性分析由以上时程图可以得到当0.1, 0.3, 0.4时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；当0.41,0.42, 0.45时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长的增加，计算结果发散得越来越快。由稳定条件，而rad/s,，由稳定条件知，当0.4时结果应当是稳定的，而且是发散与收敛的临界点，即0.4为临界步长。从图中可以看出，此时的计算结果虽然没有发散，但是结果有明显的波动，这预示着若再增大步长，则计算结果会发散。以上计算分析，说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。专心-专注-专业