I am finding response of multi degree of freedom system by using Newmarks Beta method. But i am getting some unrealistic value of response.Follwoing is my code.Please help me . Thank you in advance.
u=zeros(24,1); % assume displacemnt
du=zeros(24,1); % assume velocity
M_m=transpose(V_linear)*M*(V_linear); % Modal mass
C_m=transpose(V_linear)*C*(V_linear); % Modal damping
K_m=transpose(V_linear)*K_linear_1*(V_linear); % Modal stiffaness
q=(transpose(V_linear)*M*(V_linear))\transpose(V_linear)*M*u; % Modal displacement ,A\B == inv(A)*B
dq=(transpose(V_linear)*M*(V_linear))\transpose(V_linear)*M*du; % Modal velocity
P(:,1)=transpose(V_linear)*F(:,1); % Modal force
ddq=(M_m)\(P(:,1) - C_m*dq(:,1) - K_m*q(:,1)); % Modal acceleration
gamma=1/2;
beta=1/4;
K_cap=K_linear_1+ ((gamma)/(beta*dt))*C_m+(1/(beta*dt^2))*M_m;
a_n=(1/(beta*dt))*M_m +(gamma)/(beta)*C_m;
b=(1)/(2*beta)*M_m+dt*(gamma/(2*beta)-1)*C_m;
for i=1:(length(time)-1)
P(:,i+1)=transpose(V_linear)*F(:,i+1);% Modal force
P_cap=P(:,i+1)-P(:,i)+(a_n)*dq(:,i)+b*ddq(:,i);
delq=(K_cap)\P_cap; % A\B == inv(A)*B
deldq=(gamma)/(beta*dt)*delq-(gamma/beta)*dq(:,i)+dt*(1-gamma/(2*beta))*ddq(:,i);
delddq=1/(beta*dt^2)*delq -1/(beta*dt)*ddq(:,i)-1/(2*beta)*ddq(:,i);
q(:,i+1)=q(:,i)+delq;
dq(:,i+1)=dq(:,i)+deldq;
ddq(:,i+1)=ddq(:,i)+delddq;
u(:,i+1)=(V_linear)*q(:,i+1);
du(:,i+1)=(V_linear)*dq(:,i+1);
ddu(:,i+1)=(V_linear)*ddq(:,i+1);
tu=u(:,i+1);
tq=q(:,i+1);
tdq=dq(:,i+1);
tddq=ddq(:,i+1);
error=1;
[K_nonlinear]= K_nonlinear_coeeficient_calulcation(u,K_linear_1);
while (error>=0.001)
ttu=tu;
ttq=tq;
ttdq=tdq;
ttddq=tddq;
Pti(:,i+1)=transpose(V_linear)*F(:,i+1) - transpose(V_linear)*(K_nonlinear)*(V_linear)*ttq;
Pcapi=Pti(:,i+1)-P(:,i)+(a_n)*ttdq + b*ttddq;
delqi=K_cap\Pcapi;
deldqi=(gamma/(beta*dt))*delqi-(gamma/beta)*dq(:,i)+dt*(1-gamma/(2*beta))*ddq(:,i);
delddqi=1/(beta*dt^2)*delqi -1/(beta*dt)*ddq(:,i)-1/(2*beta)*ddq(:,i);
tq=q(:,i)+delqi;
tdq=dq(:,i)+deldqi;
tddq=ddq(:,i)+delddqi;
tu=(V_linear)*tq;
error=((tu-ttu)/ttu)*100;
end
u(:,i+1)=tu;
du(:,i+1)=(V_linear)*tdq;
ddu(:,i+1)=(V_linear)*tddq;
end
