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

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

2、种组合来表示，将微分方程组的求解问题转化为代 数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反响。中心差分法是一种显示的积分法，它基于用有限差分代昔位移对时间的求导(即速度和加速度)。如果 采用等时间步长，At, = At (At为常数)，那么速度与加速度的中心差分近似为=巴+.妇(1)12At 廿=Ui+i-2U|+Ui-i用u表示位移，离散时间点的运动为：Uf =u(q), iif = Uj = u(tt) (i = 0,1,2.)体系的运动方程为inii(t) + cu(t) 4-ku(t) = P(t)(3)将速度和加速度的差分近似公式和代入(3)

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

4、=U(0)l(5)uo=n(O)J3根据初始条件来确定根据中心差分公式O一仲.U1-2U0+U-1w& = 消去U1得到U-1的公式：u_i =UQ Atu0+ u0结构动力学中心差分法1其中零时刻加速度值论可以由 t = 0 时的运动方程得到即蜘=0一两。一如。)这样就町以根据初始条件得到 UT，然后再将初始条件应用于公式(4)中，逐步求出不同时刻的运动。中心差分法分析时的具体计算步骤：(1)根本数据准备与初始条件计算己知：初始位移口 0、&和初始荷载值 Po来计算廿 0和足. 1/一 1=UQAtUo 4- -TUQ(2)计算等效刚度和中心差分法计算公式中的系数- m c

5、k =z-dAt22At2m a = kAt2m cb =-At22At因此中心差分法计算公式可以表示为：心=R - a% -bin(3)根据4及以前的时刻的运动求得4+1时刻的运动Pf=Pt- aut_蚓_巧+i = 2/k(4)下一步计算中用i+1代替i,对于线弹性体系重复第3步计算步骤，对于非线性弹性体系，重复第2和第3计算步骤。以上的中心差分法逐步计算公式具有2阶精度，即误差ocO(At2)；并且是有条件稳定的，稳定条件为：At 7T式中，7；为结构的自振周期，对于多自由度体系那么为结构的最小自振周期。算例本算例根据结构动力学 48 页算例 3.1 数据编写，稳定条件为dt = 0.1

6、6s对于一个单层框架结构，假设楼板刚度无限大，且结构质景:集中于楼层，其质= 9240kg.刚度K=1460KN/m、阻尼系数C = 6.41KN s/m,对结构施加动力荷载P= 73000 siiiO.5TTt假设结构处于线弹性 状态，用中心差分法计算结构的自由振动反响。采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u0=0.05m和初速度v0 = 0取不同T的步长分别计算，以弗证中心差分法的稳定条件也o(8)T 7先计算由稳定条件At -=,而/结构动力学中心是分法23当dt = 0.2s时，可以得到如下提示:不满足检定条件dt=Tn/pi,清重新输入符合稳定条件的时间步

7、长dtoT2那么函 c=A0(1.3).u0=A0(1.4);v0=A0(1.5).all_time=A0(1.6);P0=A0(l7),dt=A0(l8);ifdt2*sqrt(nVk)%判断时间步长是否满足稳定条件disp(不满足稳定条件 dtcEn/pi,清重新输入符合税定条件的时间步 K 此returnelsaf 0dt=2asqrt(m/k)disp(涡足检定条件为：dtv=zeros(size(t) ac=zeros(size(t); u(:,2)=u0,v(:,2)=v0,ac(:,2)=(P0-c*v(:.2)-k*u(:.2)/m;u(:,l)=u(:,2)-dt*v(,2)

8、+(dt)A2)*ac(,2)/2.ek=m/(dt)+c/(2 dt); a=k.(2 m)/(dt)，b=m/(dt)c/(2Mt);p(:,2)=PO*sin(0);ep(:,2)=pC.2)-a*uC.2)-b*uC.l),u(:,3)=ep(:,2)/ek, for i=3 :nnp(:,i)=PO*sin(5*pi*(i-2)*dt),q)(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-l); %得出所需要结果u(:,i+l)=ep(:,i)/ek,v(:,i)=(u(:j+l)-u(ti-l)y(2*dt).endt=t(:,lend-l),u=u(,2:end-l)

9、.v=v(:,2:end),ac=ac(:,2:enci),p=p(t2 end);ep=ep(t2 end),%. 绘制位移、速度、加速度时程曲线.%plot(tu.*b-oO.hold oaplotv/g-pXholdon.plotacjr x),gnd on,xlabel(fuj(s)*)ylabel(,位移(m)速度(m/s)加速度(m/sA2),)title(,顶层 u,v,ac 的时程曲线)，subplot(3JJ),plot(tu,l)-9，grid,xlabelCRj(s)9.ylabel(i(m)9，tide(,(i 的时程曲线),legend(位移u1)subplot(3J

10、J).Flot(tv.V).gnd.xlabel(,Hr 间(sylabel(魅度(m/s)，title(，速度 v 的时程曲线legend(速度 V) subplot。,1.3),ploKLacVXgnKxlabeK时间(s),ylabel(加速度加速度ac 的时程曲线 OJcgendf加速度 ac1)%将时间分步，采用等时间步长：%计算 t 的向魅长度，得出步数：%设定存储 U的矩阵：%设定存储 V的矩阵：%设定存储 ac 的矩阵：%赋值向量第 2项为讪；%赋值向虽第 2项为 vO：%求出初始加速度 acO：%计算初始条件项：%计算等效刚度：%计算方程系数：%给出初始荷载条件：%计算初始等

11、效荷载：%计算位移 ul=u(,3)%从第二项开始进行中心是分法计算%给出荷载条件，按照箭谐荷裁计尊：%计尊位移虽:%计算速度技；ac(ti)=(u( .1+1 )-2su( .i)+u(: .1-1 )/(dtA2)t% 计算加速度结构动力学中心差分法4中心差分法求解单自由度体系的自由振动问题算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质帛M=2000kg、刚度K=50KN/m、阻尼系数C=3KNs/m,假设结构处于线弹性状态，用中心差分法计算结构的自由振动反响。采用MATLAB语言编程，并以单自由度体系为例进行计算，设初位移uO = O和初速度vO=O,取不同的步

12、长分T别计算，以验证中心差分法的稳定条件Atn先计算At,由稳定条件At =,而% = JE0000= 5函人,n a)nv M v 2000T ，)那么At = = - = 0.4兀牡5所以本次计算取也=0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算MATLAB 程序清单function u/v/ac=centraldifferent(M/C/K/uO/vO/time/dt)%本程序采用中心差分法计算结构的动力响应%本程序是既町以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态:%-% 输入参数%-% M-质量短阵% c-阻尼矩阵% K-刚度

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

14、C/2*dt% Pe=-a*X(t)-b*X(t-dt)结构动力学中心是分法5% 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)/(dtA2)%-初始条件-% XO,=(-C*XO-K*XO)% X(-l)=X0-X0,*dt+X0*(dtA2)/2%.Copyright by zhouhuaping(S202104232)clear allM=inputf输入质量矩阵M : C司nputf输入阻尼矩阵C：)； K=input(输入刚度矩阵K： )； uO=inputC输入初始位移uO：vO=in

15、putf输入初始速度vO：); time=input(输入模拟时间time：); dt=input(输入时间步长dt:)；m,m=size(K);n=time/dt;u=zeros(mzfloor(n)+l);v=zeros(m,floor(n)+l); ac=zeros(m/floor(n)+l); P=zeros(m,floor(n)+1); u(:,2)=u0; v(:,2)=v0;Ke=M/(dt)+(C)/(2*dt);a=K-2*M/dtA2;b=M/dtA2-C/(2*dt);for i=3:l rfloor(n)+l;t=(i-2)*dt;ac(:/2)=M(-K*u(:/2)

16、-C*v(:/2)；u(:,l)=u(:,2)-v(:,2)*dt+(ac(:,2)*(dt)/2;Pe= -a * u(:,i-l) -b* u(:,i-2);u(:J)=KePe;v(:J)=(u(:J)-u(:J-2)/(2*dt);ac(:J)=(u(:,i)-2*u(:J-l)+u(:J-2)/(dtA2);end%- %绘制位移、速度、加速度时程曲线%-t=0:dt:time;subplot(2,2,l),plot(t,u(m,:),k-),grid,xlabel(时间(s),yiabel(位移(顶层位移的时程曲线subplot(2,2,2),plot(t,v(m,:),T-),g

17、rid,xlabek时间),ylabel(速度(m/s),title(顶层速度的时程曲线)；subplot(2,2,3),plot(t,ac(m,:)b。grid,xlab时间(s)ylabek加速度(m/sA2)*),titieC顶层加速度的时程曲 线。；%- end运行centra Idifferent.M文件 输入参数：K=50000; M=2000; C=3000; uO=O; v0=0;time = 20s： dt= ?中心差分法计算结果稳定性分析由以上时程图可以得到当&=0.1, 0.3, 0.4时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳 定的；当At=0.41,0.42, 0.45时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长t的增加，计算结果发散得越来越快。由稳定条件At =,而= 区=WW=5rad/s,兀conVM V 2000T 2%计算步数%设定存储位移矩阵%设定