%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%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
