4NewMark linear acceleration.zip_Newmark 位移_acceleration_newmark

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%clear %%%clc %% newmark method linear acceleration method Gamma = 1/2; Beta = 1/6; %% constant number M = 0.2533; K = 10; C = 0.1592; %% initial calculations MM=1; [HC,dt]=K_L(MM); %%%HC = 0.1; % hardening coefficient %%%dt = 0.1; % step lenth t = [0.0:dt:0.6]; p = 10*sin(pi*t/0.6); t_add = [0.6+dt:dt:10]; T = [t,t_add]; p = [p,t_add*0]; N = size(T,2); U = zeros(1,N); V = zeros(1,N); A = zeros(1,N); fs = zeros(1,N); Kt = zeros(1,N); PP = zeros(1,N); fx = zeros(1,N); fy = zeros(1,N); fm = zeros(1,N); fr = zeros(1,N); fq = zeros(1,N); ft = zeros(1,N); %% first calculation c1=1/(Beta*dt^2); c2=Gamma/(Beta*dt); c3=1/(Beta*dt); c4=1/(2*Beta)-1; c5=Gamma/Beta-1; c6=(Gamma/(2*Beta)-1)*dt; c7=(1-Gamma)*dt; c8=Gamma*dt; a1=M*c1+c2*C; a2=M*c3+c5*C; a3=c4*M+c6*C; U(1) = 0; fs(1) = 0; Kt(1) = 10; fy(1) = 7.5; fm(1) = 0; ut = 0; ux = 0; fr(1) = 0; ur = 0; uq = 0; fq(1) = 0; up = 0; %% lterative calculation for i = 1:N-1 U(i+1) = U(i); fs(i+1) = fs(i); Kt(i+1) = Kt(i); pp(i+1) = p(i+1)+a1*U(i)+a2*V(i)+a3*A(i); %% Internal circulation while(true) R = pp(i+1)-fs(i+1)-a1*U(i+1); if abs(R)<1e-3 break; end KKt = Kt(i+1)+a1; du = KKt\R; U(i+1) = U(i+1) + du; fs(i+1) = fs(i) + Kt(i) * (U(i+1)-U(i)); %% State of determination if fm(i) == 0 fm(i) = fy(1); fr(i+1) = fr(i); end if (fs(i+1) < fm(i)) && (fs(i+1) >= 0) Kt(i+1) = Kt(1); fm(i+1) = fm(i); fr(i+1) = fr(i); ft(i+1) = ft(i); %%%%第二段 elseif (fs(i+1) >= fm(i) && fs(i+1) >0) && (U(i+1)-U(i)) >= 0 if ut == 0; ut = U(i); end Kt(i+1) = Kt(1)*HC; fs(i+1) = fm(i)+(U(i+1)-ut)*Kt(1)*HC; fm(i+1) = fm(i); fr(i+1) = fr(i); %%%%第三段 elseif (fs(i+1) >= fm(i) && fs(i+1) >0) && (U(i+1)-U(i)) <= 0 if fr(i) == 0 && up == 0 fr(i) = fs(i+1); up = U(i); end Kt(i+1) = Kt(1); fr(i+1) = fr(i); fs(i+1) = fr(i) + Kt(1) * (U(i+1)-up); %%%%第四段 elseif fs(i+1) > -(2*fy(1)-fr(i)) && ((U(i+1)-U(i)) < 0) Kt(i+1) = Kt(1); fr(i+1) = fr(i); fm(i+1) = fm(i); fs(i+1) = fs(i) + Kt(1) * (U(i+1)-U(i)); %%%%第五段 elseif fs(i+1) <= -(2*fy(1)-fr(i)) && (U(i+1)-U(i)) <= 0 if ux == 0 ux = U(i); end Kt(i+1) = Kt(1)*HC; fs(i+1) = -(2*fy(1)-fr(i))+(U(i+1)-ux)*Kt(1)*HC; fr(i+1) = fr(i); fm(i+1) = fm(i); %%%%第六段 elseif fs(i+1) <= -(2*fy(1)-fr(i)) && ((U(i+1)-U(i)) > 0) if uq == 0 uq = U(i); fq(i) = fs(i+1); end Kt(i+1) = Kt(1); fs(i+1) = fq(i) + (U(i+1)-uq) * Kt(1); fq(i+1) = fq(i); ft(i) = 2*fy(1)+fq(i); ft(i+1) = ft(i); fm(i+1) = ft(i); fr(i+1) = fr(i); elseif fs(i+1)<0 && fs(i+1)-fs(i)>0 fs(i+1) = fs(i) + Kt(1) * (U(i+1)-U(i)); Kt(i+1) = Kt(1); fm(i+1) = ft(i); fr(i+1) = fr(i); end end %% Calculating speed and acceleration V(i+1) = c2*(U(i+1)-U(i))-c5*V(i)-c6*A(i); A(i+1) = c1*(U(i+1)-U(i))-c3*V(i)-c4*A(i); end %% end Calculation %% plot figure figure(1) plot(T,U(1,:),'g-.','LineWidth',2) xlabel('times:s') ylabel('Displacement:ft') legend('NewMark linear') hold on grid on figure(2) plot(T,V(1,:),'g-.','LineWidth',2) xlabel('times:s') ylabel('Velocity:ft/s') legend('NewMark linear') hold on grid on figure(3) plot(T,A(1,:),'g-.','LineWidth',2) xlabel('times:s') ylabel('Acceleration:ft/s^2') legend('NewMark linear') hold on Z=zeros(4,N); Z=[U;fs;V;A]; grid on figure(4) plot(U(1,:),fs(1,:),'g-.','LineWidth',2) legend('NewMark linear'); hold on