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

中心差分法计算单自由度体系动力反映的报告前言基于叠加原理的时域积分法与频域积分法一样，都假设结构在在全部反应过程中都是线性的。而时域逐步积分法只是假设结构本构关系在一个微小的时间步距内是线性的，相当于分段直线来逼近实际的曲线。时域逐步积分法是结构动力问题中研究并应用广泛的课题。中心差分法是一种目前发展的一系列结构动力反应分析的时域逐步积分法的一种，时域逐步积分法还包括分段解析法、平均常加速度法、线性加速度法、Newmarket−β和Wilson−θ法等。中心差分法(centraldifferencemethod)原理中心差分法的基本思路将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。中心差分法是一种显示的积分法，它基于用有限差分代替位移对时间的求导（即速度和加速度）。如果采用等时间步长，∆t𝑖=∆t（∆t为常数），则速度与加速度的中心差分近似为u̇𝑖=𝑢𝑖+1+𝑢𝑖−12∆𝑡(1)𝑢̈𝑖=𝑢𝑖+1−2𝑢𝑖+𝑢𝑖−1∆𝑡2(2)用u表示位移，离散时间点的运动为：u𝑖=𝑢(𝑡𝑖)，u̇𝑖=𝑢̇(𝑡𝑖)，ü𝑖=𝑢̈(𝑡𝑖)(𝑖=0,1,2…)体系的运动方程为m𝑢̈(𝑡)+𝑐𝑢̇(𝑡)+𝑘𝑢(𝑡)=𝑃(𝑡)(3)将速度和加速度的差分近似公式(1)和(2)代入(3)中得出在𝑡𝑖时刻的运动方程，将方程整理得到𝑢𝑖+1由𝑢𝑖和𝑢𝑖−1表示的两步法的运动方程(4)：(𝑚∆𝑡2+𝑐2∆𝑡)𝑢𝑖+1=𝑃𝑖−(𝑘−2𝑚∆𝑡2)𝑢𝑖−(𝑚∆𝑡2−𝑐2∆𝑡)𝑢𝑖−1(4)这样就可以根据𝑡𝑖及以前的时刻的运动求得𝑡𝑖+1时刻的运动。中心差分法属于两步法，用两步法计算时存在起步问题，必须要给出相邻的两个时刻的位移值，才能逐步计算。对于地震作用下结构的反应问题和一般的零初始条件下的动力问题，可以用(4)直接计算，因为总可以假设初始的两个时间点（一般取𝑖=0,−1）的位移等于零。但是对应于非零初始条件或零时刻外荷载很大时，需要进行一定的分析，建立两个起步时刻（即𝑖=0,−1）的位移值。假设给定的初始条件为𝑢0=𝑢(0)𝑢̇0=𝑢̇(0)}(5)根据初始条件来确定𝑢−1。根据中心差分公式𝑢̇0=𝑢1+𝑢−12∆𝑡𝑢̈0=𝑢1−2𝑢0+𝑢−1∆𝑡2}(6)消去𝑢1得到𝑢−1的公式：𝑢−1=𝑢0−∆𝑡𝑢̇0+∆𝑡22𝑢̈0(7)结构动力学中心差分法1其中零时刻加速度值𝑢̈0可以由t=0时的运动方程得到即𝑢̈0=1𝑚(𝑃0−𝑐𝑢̇0−𝑘𝑢0)(8)这样就可以根据初始条件得到𝑢−1，然后再将初始条件应用于公式(4)中，逐步求出不同时刻的运动。中心差分法分析时的具体计算步骤：(1)基本数据准备与初始条件计算已知：初始位移𝑢0、𝑢̇0和初始荷载值𝑃0来计算𝑢̈0和𝑢−1𝑢̈0=1𝑚(𝑃0−𝑐𝑢̇0−𝑘𝑢0)𝑢−1=𝑢0−∆𝑡𝑢̇0+∆𝑡22𝑢̈0(2)计算等效刚度和中心差分法计算公式中的系数𝑘̂=𝑚∆𝑡2+𝑐2∆𝑡𝑎=𝑘−2𝑚∆𝑡2b=𝑚∆𝑡2−𝑐2∆𝑡因此中心差分法计算公式可以表示为：𝑘̂𝑢𝑖+1=𝑃𝑖−𝑎𝑢𝑖−𝑏𝑢𝑖−1(3)根据𝑡𝑖及以前的时刻的运动求得𝑡𝑖+1时刻的运动𝑃̂𝑖=𝑃𝑖−𝑎𝑢𝑖−𝑏𝑢𝑖−1𝑢𝑖+1=𝑃̂𝑖𝑘̂⁄(4)下一步计算中用𝑖+1代替𝑖，对于线弹性体系重复第3步计算步骤，对于非线性弹性体系，重复第2和第3计算步骤。以上的中心差分法逐步计算公式具有2阶精度，即误差ϵ∝O(∆𝑡2)；并且是有条件稳定的，稳定条件为：∆𝑡≤𝑇𝑛𝜋式中，𝑇𝑛为结构的自振周期，对于多自由度体系则为结构的最小自振周期。算例本算例根据结构动力学48页算例3.1数据编写，稳定条件为𝐝𝐭=𝟎.𝟏𝟔𝐬对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=9240kg、刚度K＝1460KN/m、阻尼系数C＝6.41KN∙s/m，对结构施加动力荷载P=73000∙sin0.5𝜋𝑡假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u0＝0.05m和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件nTt。先计算t，由稳定条件nnTt2，而146000012.579240nKMrad/s，结构动力学中心差分法2则20.16nnTts，所以本次计算取t＝0.2，0.1，0.05分别进行计算。计算结果与分析1）当dt=0.05s时，可以得到位移u，速度v，加速度ac的时程曲线如下：2）当dt=01s时，可以得到位移u，速度v，加速度ac的时程曲线如下：3）当dt=0.2s时，可以得到如下提示：不满足稳定条件:dt=Tn/pi,请重新输入符合稳定条件的时间步长dt。结构动力学中心差分法3附录%m=质量；k=刚度；c=阻尼；u0=初始位移；v0=初始速度；all_time=所用时间；P0=荷载幅值；dt=时间步长；%u=位移；v=速度；ac=加速度；ek=等效刚度；p=荷载；ep=等效荷载；t=时间；clear%A0=input('请按格式和顺序输入初始矩阵，如A0=[m,k,c,u0,v0,all_time,P0,dt]=[9240146000064100.050307300*30.05]=');A0=[9240146000064100.05020730000.05];m=A0(1,1);k=A0(1,2);c=A0(1,3);u0=A0(1,4);v0=A0(1,5);all_time=A0(1,6);P0=A0(1,7);dt=A0(1,8);ifdt2*sqrt(m/k)%判断时间步长是否满足稳定条件disp('不满足稳定条件:dt=Tn/pi,请重新输入符合稳定条件的时间步长dt')returnelseif0dt=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(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/(dt^2)+c/(2*dt);%计算等效刚度；a=k-(2*m)/(dt^2);b=m/(dt^2)-c/(2*dt);%计算方程系数；p(:,2)=P0*sin(0);%给出初始荷载条件；ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1);%计算初始等效荷载；u(:,3)=ep(:,2)/ek;%计算位移u1=u(:,3)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(:,i)=(u(:,i+1)-u(:,i-1))/(2*dt);%计算速度量；ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1))/(dt^2);%计算加速度量；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'),holdon,plot(t,v,'g--p'),holdon,plot(t,ac,'r:x'),gridon,xlabel('时间(s)'),ylabel('位移(m)速度(m/s)加速度(m/s^2)'),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的时程曲线');legend('速度v')subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),title('加速度ac的时程曲线');legend('加速度ac')结构动力学中心差分法4中心差分法求解单自由度体系的自由振动问题算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量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----------------