ممنون- در خط 61 for fybar=0.1:0.1:1 اصلا کار نمیکندو fybar دایما با مقدار 0.1 کار میکند
کد:
P=input('please input data matrix: ');
ksi=input('please input damping ratio: ');
m=input('please input the mass of the system: ');
k=input('please input the stiffness of the system: ');
p=P*(-m);
u0=0;
udot=0;
udob=0;
A0=0;
G=0.5;
B=0.25;
t=0.02;
count = 1;
for Tn=0.01:0.2:50
c=2*m*ksi*2*pi/Tn;
a1 = (m/(B*(t^2))) + ((G*c)/(B*t));
a2 = (m/(B*t)) + (( (G/B)-1) *c);
a3 = m * ((1/(2*B))-1) + (c * t * ((G/(2*B))-1));
K = k+a1;
z = size(P)-2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Uo %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
udob(1,1)=0;
udot(1,1)=0;
u(1,1)=0;
for i=1:z
pd(i+1,1)=p(i+1,1)+(a1*u(i,1))+(a2*udot(i,1))+(a3*udob(i,1));
u(i+1,1)=pd(i+1,1)/K;
udot(i+1,1)=((2/t)*(u(i+1,1)-u(i,1)))-udot(i,1);
udob(i+1,1)=(4/(t^2)) * (u(i+1,1)- u(i))- ((4/t)*udot(i+1,1)) - udob(i,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Calculation of Um %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Uo=max(u);
fo=k*Uo;
disp(fo);
kd=k;
delu=0;
%nonlinear analysis( New mark method)
%iterations used
%first iteration n=1
%h=i
%
final(count,1)=Tn;
final(count,2)=Uo;
countfy=3;
for fybar=0.1:0.1:1
fy=fybar*fo;
uy=fy/k;
%% section Um
for h=1:z
for n=1:100
fs(h+1,1)=(m*udob(h+1,1))+(c*udot(h+1,1))+(p(h+1,1));
R(h+1,1)=pd(h+1,1)-fs(h+1,1)-(a1*u(h+1,1));
if R(h+1,1)<0.001
un(h+1,1)=fs(h+1,1)/k;
else
delu=R(h+1,1)/kd;
un(h+1,1)=u(h+1,1)+delu;
fsn(h+1,1)=fs(h,1)+(k*(un(h+1,1)-un(h,1)));
if fsn(h+1,1)>fy
fs(h+1,1)=fy;
un(h+1,1)=fs(h+1,1)/k;
kd = a1;
else
un(h+1,1)=fsn(h+1,1)/k;
kd = kd + a1;
end
end
end
end
Um = max(un);
final(count,countfy)=Um;
countfy = countfy + 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
count = count + 1;
end
%% moo calculation
for i=1:250
fyyy=0.01;
for j=1:10
moo(i,j) = final(i,j+2)/(final(i,2)*fyyy);
fyyy = fyyy + 0.1;
end
end
%%
MooValue = input('Enter ductility factor: ');
n1 = 0;
n2 = 0;
fybarM=0.1
result =0;
for i=1:250
for j=1:9
n1=moo(i,j);
n2=moo(i,j+1);
if (MooValue <= n2)
if (MooValue > n1)
fybarM(i,1) = j * (0.1);
end
end
end
end
count = 1;
for i=0.01:0.2:50
result(count, 1) = i; % Tn
result(count, 2) = fybarM(count,1); % Fy bar
result(count, 3) = fybarM(count,1)*final(count, 2); % Uy
result(count, 4) = result(count,3)*2*pi/i; % Vy
result(count, 5) = result(count,3)*((2*pi/i)^2); % Ay
count = count + 1;
end
%