《中心差分法在单自由度中的应用.doc》由会员分享，可在线阅读，更多相关《中心差分法在单自由度中的应用.doc（5页珍藏版）》请在金锄头文库上搜索。
1、北京工业大学建筑工程学院 结构动力学中心差分法求解单自由度体系的自由振动问题前言时域逐步积分法是根据运动方程，引进某些假设，建立由t时刻状态向量、到t时刻的状态向量、的递推关系，从而从t0时刻的初始状态向量、出发，逐步求出各时刻的状态向量，由于引进的假设条件不同，可以有各种不同的方法，下面主要介绍一种时域逐步积分方法中心差分法。中心差分法(central difference method)原理1 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反
2、应。中心差分法只在相隔一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，则速度与加速度的中心差分近似为： (a) (b)而离散时间点的运动为 ( 0,1,2,3,)由体系运动方程为： (c)将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到时刻的运动方程： （d）在（d）式中，假设和是已知的，即在及以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到： (e)由式（e）就可以根据及以前时刻的运动，求得时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。假设给定的初始条件为 （g）由式（g）确定。在零时
3、刻速度和加速度的中心差分公式为： （h） （i）将式（i）消去得： （j） 而零时刻的加速度值可以用t0时的运动方程 确定即 （k）这样就可以根据初始条件和初始荷载，就可以根据上式确定的值。下面给出采用中心差分法分析时的具体计算步骤：（1） 基本数据准备和初始条件计算 （2） 计算等效刚度和中心差分计算公式中的相关系数 （3） 根据及以前时刻的运动，计算时刻的运动 如果需要，可计算（4）下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。以上为中心差分法逐步计算公式，其具有2阶精度，即误差；并且为有条件稳定，稳定条件为： 上式中，为结构的自振周期，
4、对于多自由度结构体系则为结构的最小自振周期。算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg、刚度K50KN/m、阻尼系数C3KNs/m，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。采用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,
5、u0,v0,time,dt) % 本程序采用中心差分法计算结构的动力响应% 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态；% -%输入参数%-% M-质量矩阵% 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
6、(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-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 zhouhuapi
7、ng(S201004232)-clear allM=input(输入质量矩阵M ：);C=input(输入阻尼矩阵C：);K=input(输入刚度矩阵K：);u0=input(输入初始位移 u0： );v0=input(输入初始速度 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=
8、zeros(m,floor(n)+1); %设定存储荷载矩阵 u(:,2)=u0; %给定初位移v(:,2)=v0; %给定初速度Ke=M/(dt2)+(C)/(2*dt); %等效刚度Ke及系数a、ba=K-2*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(:,
9、i)=KePe; %计算位移 v(:,i)=(u(:,i)-u(:,i-2)/(2*dt); %计算速度 ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2)/(dt2); %计算加速度end%-%绘制位移、速度、加速度时程曲线%- 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(顶层速度的时程
10、曲线); subplot(2,2,3),plot(t,ac(m,:),b-),grid,xlabel(时间(s),ylabel(加速度(m/s2),title(顶层加速度的时程曲线); %-end运行centraldifferent.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为临界步长。从图中可以看出，此时的计算结果虽然没有发散，但是结果有明显的波动，这预示着若再增大步长，则计算结果会发散。以上计算分析，说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。第5页
