Hello,
I had previously asked a question about my code. I have corrected the formulas but I still have a problem. My problem about this code is "a1" values. They are calculated correctly, except for the ones computed for x=Lb+d, 2Lb+d, ... 6Lb+d, which are limits of some if statements. For example, the 23th value of a1 (corresponds to x=Lb+d) should've been 12.35 but instead, somehow it is calculated zero. b1 show have been zero but instead, it is 12.35. I can't see what I am doing wrong. Any help is appriciated.
% MR-IL & Response
sp=7; % no of spans
d1=1.95; % distance between the unit load and point x
x1=zeros(sp*N+1,1); d_xs2=zeros(1,sp*N+1);
b1=zeros(sp*N+1,1); a1=zeros(sp*N+1,1); k1=226833; k2=226833; % [calculated-kN/m]
AAB=zeros(1,sp*N+1); % Area multiplication of AB triangle
ABC=zeros(1,sp*N+1); % Area multiplication of BC triangle
ACD=zeros(1,sp*N+1); % Area multiplication of CD triangle
d_xs1=zeros(1,sp*N+1); % Bending component of MR-IL
d_xs=zeros(1,sp*N+1); % Total MR-IL
resp_x=zeros(1,sp*N+1); % MR-IL Response
for i=1:sp*N+1
x1(i)=(i-1)*(Lb/N);
if x1(i)<d1 % before load arrives
b1(i)=0;
a1(i)=0;
d_xs2(1,i)=0; % support deflection component of MR-IL
AAB(1,i)=0;
ABC(1,i)=0;
ACD(1,i)=0;
d_xs1(1,i)=0;
d_xs(1,i)=0;
elseif x1(i)==d1 % when the load is on the first support
b1(i)=0;
a1(i)=0;
d_xs2(1,i)=(Lb-x1(i))/(k1.*Lb);
AAB(1,i)=0;
ABC(1,i)=0;
ACD(1,i)=0;
d_xs1(1,i)=0;
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (d1<x1(i)) && (x1(i)<Lb) % load and x in the same span
a1(i)=(i-4)*L;
b1(i)=Lb-a1(i);
if a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + (x1(i,1)*(b1(i)-a1(i))/(k1*Lb^2));
elseif a1(i) > b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + (x1(i,1)*(a1(i)-b1(i))/(k1*Lb^2));
end
AAB(1,i)=(b1(i)*(b1(i)-d1)*a1(i)^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (Lb<=x1(i)) && (x1(i)<=(Lb+d1)) %load is one span behind point x.
a1(i)=(i-4)*L;
b1(i)=Lb-a1(i);
d_xs2(1,i) =(2*Lb-x1(i))*a1(i)/(k1*Lb^2);
d_xs1(1,i)=0; % load is in the previous span so, no bending effect occurs in this span
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (Lb+d1<x1(i)) && (x1(i)<2*Lb)
a1(i)=(i-4)*L-Lb;
b1(i)=Lb-a1(i);
if a1(i)>b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-Lb)/(k1*Lb^2));
elseif a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (2*Lb<=x1(i)) && (x1(i)<=(2*Lb+d1))
a1(i)=(i-4)*L-Lb;
b1(i)=Lb-a1(i);
d_xs2(1,i) = (3*Lb-x1(i)).*a1(i)./(k1*Lb^2);
d_xs1(1,i)=0; % load is in the previous span so, no bending effect occurs
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (2*Lb+d1<x1(i)) && (x1(i)<3*Lb)
a1(i)=(i-4)*L-2*Lb;
b1(i)=Lb-a1(i);
if a1(i) > b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-2*Lb)/(k1*Lb^2));
elseif a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-2*Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (3*Lb<=x1(i)) && (x1(i)<=(3*Lb+d1))
a1(i)=(i-4)*L-2*Lb;
b1(i)=Lb-a1(i);
d_xs2(1,i) =(4*Lb-x1(i)).*a1(i)./(k1*Lb^2);
d_xs1(1,i)=0;
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (3*Lb+d1<x1(i)) && (x1(i)<4*Lb)
a1(i)=(i-4)*L-3*Lb;
b1(i)=Lb-a1(i);
if a1(i)>b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-3*Lb)/(k1*Lb^2));
elseif a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-3*Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (4*Lb<=x1(i)) && (x1(i)<=(4*Lb+d1))
a1(i)=(i-4)*L-3*Lb;
b1(i)=Lb-a1(i);
d_xs2(1,i) = (5*Lb-x1(i)).*a1(i)./(k1*Lb^2);
d_xs1(1,i)=0;
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (4*Lb+d1<x1(i)) && (x1(i)<5*Lb)
a1(i)=(i-4)*L-4*Lb;
b1(i)=Lb-a1(i);
if a1(i)>b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-4*Lb)/(k1*Lb^2));
elseif a1(i)< b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-4*Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (5*Lb<=x1(i)) && (x1(i)<=(5*Lb+d1))
a1(i)=(i-4)*L-4*Lb;
b1(i)=Lb-a1(i);
d_xs2(1,i) = (6*Lb-x1(i)).*a1(i)./(k1*Lb^2);
d_xs1(1,i)=0;
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (5*Lb+d1<x1(i)) && (x1(i)<6*Lb)
a1(i)=(i-4)*L-5*Lb;
b1(i)=Lb-a1(i);
if a1(i)>b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-5*Lb)/(k1*Lb^2));
elseif a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-5*Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
elseif (6*Lb<=x1(i)) && (x1(i)<=(6*Lb+d1))
a1(i)=(i-4)*L-5*Lb;
b1(i)=Lb-a1(i);
d_xs2(1,i) =(7*Lb-x1(i)).*a1(i)./(k1*Lb^2);
d_xs1(1,i)=0;
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i));
elseif (6*Lb+d1<x1(i)) && (x1(i)<=7*Lb)
a1(i)=(i-4)*L-6*Lb;
b1(i)=Lb-a1(i);
if a1(i)>b1(i)
d_xs2(1,i)=(b1(i)/(k1*Lb)) + ((a1(i)-b1(i))*(x1(i)-6*Lb)/(k1*Lb^2));
elseif a1(i)<b1(i)
d_xs2(1,i)=(a1(i)/(k1*Lb)) + ((b1(i)-a1(i))*(x1(i)-6*Lb)/(k1*Lb^2));
end
AAB(1,i)=(b1(i).*(b1(i)-d1).*a1(i).^3)/(3*Lb^2);
ABC(1,i)=d1*(2*a1(i).^2.*b1(i).*(b1(i)-d1) + a1(i).*b1(i).*(a1(i)+d1).*(b1(i)-d1)+(b1(i)-d1)^2.*a1(i).^2+2.*(b1(i)-d1).*2.*a1(i).*(a1(i)+d1))/(6*Lb^2);
ACD(1,i)=(a1(i).*(a1(i)+d1).*(b1(i)-d1).^3)/(3*Lb^2);
d_xs1(1,i) =(AAB(1,i)+ABC(1,i)+ACD(1,i))./(EI/1000);
d_xs(1,i) =-(d_xs1(1,i)+d_xs2(1,i)); % total MR-IL
else
disp('Error')
end
resp_x(i)= (Gvw/100).*d_xs(1,i); % response due to 1 axle
end

댓글 수: 4

Cris LaPierre
Cris LaPierre 2022년 2월 27일
"calculated dx_t values (resultant graph) are wrong for some intervals."
Which intervals?
Torsten
Torsten 2022년 2월 27일
Maybe for both if-clauses you should include an
else
disp('Error')
end
at the end in order to ensure that one of the cases has been used by MATLAB.
Sinem Tola
Sinem Tola 2022년 2월 27일
At least these ones:
(3*Lb+d1<=x1(i)) && (x1(i)<4*Lb)
(4*Lb+d1<=x1(i)) && (x1(i)<5*Lb)
(5*Lb+d1<=x1(i)) && (x1(i)<6*Lb)
(6*Lb+d1<=x1(i)) && (x1(i)<7*Lb)
Looks like the values jump up in these intervals. I cheked with an another program and results are close to the values one interval ahead of these above. For example:
(3*Lb<=x1(i)) && (x1(i)<3*Lb+d1)
(4*Lb<=x1(i)) && (x1(i)<4*Lb+d1)
(5*Lb<=x1(i)) && (x1(i)<5*Lb+d1)
(6*Lb<=x1(i)) && (x1(i)<6*Lb+d1)
Torsten
Torsten 2022년 2월 27일
@Sinem Tola comment moved here:
Thank you Torsten, I did it and it did not result with the error. Maybe I should check my formulas one more time. Because the graph should be a concave curve.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

도움말 센터File Exchange에서 Programming에 대해 자세히 알아보기

질문:

2022년 2월 27일

편집:

2022년 4월 29일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by