Newmark法-自编,newmark法原理,matlab

clc;clear %Newmark-beta法，振型分解法，振型叠加法 l=25; E=2.89e9; I=2.9; m=2303; k=1595; V=100/3.6; M=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5750]; M=diag(M); w=[]; for j=1:20%20阶 % w(j)=j^2*pi*pi*sqrt((E*I)/(m*l^4)); w(j)=(j*pi/l)^2*sqrt(E*I/m); end pk=2*k/(m*l); pf=2*5750*9.81/(m*l); nt=900;%计算步数 dt=0.001;%计算时间步长 gama=0.5; beta=0.25; a0=1/beta/dt/dt; a1=gama/beta/dt; a2=1/beta/dt; a3=1/2/beta -1; a4=gama/beta-1; a5=dt/2*(gama/beta-2); a6=dt*(1-gama); a7=dt*gama; q=zeros(21,2); v=zeros(21,2); a=zeros(21,2); tt=zeros(1,nt+1); % q=[]; % v=[]; % a=[]; k1=[]; k2=zeros(20,20); k3=[]; k4=[]; k5=[]; for i=2:(nt+1) tt(1,1)=1*dt; tt(1,i)=i*dt; t=(i-1)*dt; for i1=1:20 fai(i1)=sin((i1*pi*V*t)/l); end for b=1:20 for c=1:20 k1(b,c)=pk*fai(b)*fai(c); end end for b=1:20 k2(b,b)=w(b)*w(b); k4(b)=-k*fai(b);%车辆矩阵最后一行 k5(b,1)=-pk*fai(b);%桥梁矩阵最后一列 end k3=k1+k2; k7=[k3 k5]; k6=[k4 k]; % for b=1:21 % k7(21,b)=k6(b); % end K=[k7;k6]; f=[]; for b=1:20 f(b)=pf.*fai(b); end f=[f 0]'; Ke=K+a0.*M; Fe=f+M*(a0.*q(:,i-1)+a2.*v(:,i-1)+a3.*a(:,i-1)); q(:,i)=inv(Ke)*Fe; a(:,i)=a0.*(q(:,i)-q(:,i-1))-a2.*v(:,i-1)-a3.*a(:,i-1); v(:,i)=v(:,i-1)+a6.*a(:,i-1)+a7.*a(:,i); end u=zeros(1,nt+1); vb=zeros(1,nt+1); ab=zeros(1,nt+1); for b=1:20%振型叠加 u=q(b,:)*sin(b*pi/2)+u; vb=v(b,:)*sin(b*pi/2)+vb; ab=a(b,:)*sin(b*pi/2)+ab;%公式中x取L/2 end % u=zeros(20,nt+1); % vb=zeros(20,nt+1); % ab=zeros(20,nt+1); % for b=1:20%振型叠加 % u(b,:)=q(b,:)*sin(b*pi/2); % vb(b,:)=v(b,:)*sin(b*pi/2); % ab(b,:)=a(b,:)*sin(b*pi/2);%之前乘fai(b)，公式中x取L/2 % end % u=sum(u); % vb=sum(vb); % ab=sum(ab); qv=q(21,:); vv=v(21,:); av=a(21,:); figure; plot(tt(1,:),u) ylabel('位移/[m]'); grid on; xlabel('时间/[s]'); title('桥梁跨中位移响应时程'); figure; plot(tt(1,:),ab) ylabel('加速度/[m/s^2]'); grid on; xlabel('时间/[s]'); title('桥梁跨中加速度响应时程');