中心差分法求解单自由度体系的自由振动问题

中心差分法求解单自由度体系的自由振动问题前言时域逐步积分法是根据运动方程，引进某些假设，建立由t时刻状态向量iu、iu、iu到t＋t时刻的状态向量1iu、1iu、1iu的递推关系，从而从t＝0时刻的初始状态向量0u、0u、0u出发，逐步求出各时刻的状态向量，由于引进的假设条件不同，可以有各种不同的方法，下面主要介绍一种时域逐步积分方法－中心差分法。中心差分法(centraldifferencemethod)原理[1]中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。中心差分法只在相隔t一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，tti，则速度与加速度的中心差分近似为：tuuuii211(a)2112tuuuuiii(b)而离散时间点的运动为)(),(),(iiiiiituutuutuu(i0,1,2,3,……)由体系运动方程为：0)()()(tkutuctumi(c)将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到it时刻的运动方程：02211211iiiiiikutuuctuuum（d）在（d）式中，假设iu和1iu是已知的，即在it及it以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(iiiutctmutmkutctm(e)由式（e）就可以根据it及it以前时刻的运动，求得1it时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。假设给定的初始条件为),0(),0(00uuuu（g）由式（g）确定1u。在零时刻速度和加速度的中心差分公式为：tuuu2110（h）`210102tuuuu（i）将式（i）消去1u得：020012ututuu（j）而零时刻的加速度值0u可以用t＝0时的运动方程0000kuucum确定即)(1000kuucmu（k）这样就可以根据初始条件00,uu和初始荷载0P，就可以根据上式确定1u的值。下面给出采用中心差分法分析时的具体计算步骤：（1）基本数据准备和初始条件计算)(1000kuucmu020012ututuu（2）计算等效刚度和中心差分计算公式中的相关系数tctmk2222tmkatctmb22（3）根据it及it以前时刻的运动，计算1it时刻的运动1iibuauPkPui1如果需要，可计算tuuuii2112112tuuuuiii（4）下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。以上为中心差分法逐步计算公式，其具有2阶精度，即误差)(02t；并且为有条件稳定，稳定条件为：nTt上式中，nT为结构的自振周期，对于多自由度结构体系则为结构的最小自振周期。算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg、刚度K＝50KN/m、阻尼系数C＝3KNs/m，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u0＝0和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件nTt。先计算t，由稳定条件nnTt2，而5200050000MKnrad/s,则4.0522nnTt所以本次计算取t＝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------------质量矩阵%C------------阻尼矩阵%K------------刚度矩阵%u0-----------初始位移%v0-----------初始速度%time---------模拟时间%dt-----------时间步长%---%%%%%%输出值%%%%%%%%------%u--------------位移%v--------------速度%ac-------------加速度%-------%%%%%%%%中心差分法主要公式及原理%%%%%%%%%%-----------%MX''+CX'+KX=0%M*(X(t+dt)-2*X(t)+X(t-dt))/(dt^2)+C*(X(t+dt)-X(t-dt))/(2*dt)+K*X(t)=0%(M/dt^2+C/2*dt)*(X(t+dt))=-(K-2*M/dt^2)*X(t+dt)-(M/dt^2-C/2*dt)*X(t-dt)%-----------------等效刚度Ke等效荷载Pe和相关系数a,b-------------------------%Ke=M/dt^2+C/2*dt%a=K-2*M/dt^2%b=M/dt^2-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))/(dt^2)%------------------初始条件---------------------%X0''=(-C*X0'-K*X0)%X(-1)=X0-X0'*dt+X0''*(dt^2)/2%--------@Copyrightbyzhouhuaping(S201004232)-----clearallM=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=zeros(m,floor(n)+1);%设定存储荷载矩阵u(:,2)=u0;%给定初位移v(:,2)=v0;%给定初速度Ke=M/(dt^2)+((C)/(2*dt));%等效刚度Ke及系数a、ba=K-2*M/dt^2;b=M/dt^2-C/(2*dt);fori=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)*(dt^2))/2;%计算(0-dt)时刻位移Pe=-a*u(:,i-1)-b*u(:,i-2);%计算等效荷载Peu(:,i)=Ke\Pe;%计算位移v(:,i)=(u(:,i)-u(:,i-2))/(2*dt);%计算速度ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2))/(dt^2);%计算加速度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('顶层速度的时程曲线');subplot(2,2,3),plot(t,ac(m,:),'b-'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),title('顶层加速度的时程曲线');%----------end运行centraldifferent.M文件输入参数：K=50000;M=2000;C=3000;u0=0;v0=0;time＝20s；dt＝？中心差分法计算结果稳定性分析由以上时程图可以得到当t＝0.1,0.3,0.4时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；当t＝0.41,0.42,0.45时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长t的增加，计算结果发散得越来越快。由稳定条件nnTt2，而5200050000MKnrad/s,4.0522nnTt，由稳定条件知，当t＝0.4时结果应当是稳定的，而且是发散与收敛的临界点，即t＝0.4为临界步长。从图中可以看出，此时的计算结果虽然没有发散，但是结果有明显的波动，这预示着若再增大步长t，则计算结果会发散。以上计算分析，说明了中心差分法是有条件稳定的并验证了中心差分法的稳定条件。